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ABSTRACT 

The  aim  of  this  paper  is  to  provide  convenient  predictor- 
corrector  (P-C)  methods  for  obtaining  accurate  numerical 
solution  at  a  minimum  cost  to  first  order  ordinary  differen- 
tial equations  (ODE).   In  pursuing  this  goal,  a  unified 
development  of  the  most  popular  and  efficient  P-C  methods  is 
presented,  which  includes  derivation  of  formulas  and  analysis 
of  error  propagation  and  numerical  stability.   Each  method  is 
then  coded  and  programmed  using  the  Fortran  language.   Com- 
parative analysis  of  the  different  P-C  methods  include  both 
theoretical  and  numerical  results.   The  numerical  results 
were  obtained  by  subjecting  each  method  to  a  wide  variety  of 
test  ODE,  using  a  maximum  of  two  corrector  applications  and 
a  uniform  series  of  step  size  values.   By  systematic  compari- 
son of  the  performance  of  each  P-C  method  the  most  convenient 
P-C  sets  in  terms  of  accuracy  and  minimum  cost  are  established 
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I.   INTRODUCTION 

A  linear  first  order  ordinary  differential  equation  (ODE) 
can  be  used  as  a  mathematical  model  for  a  variety  of  phenom- 
ena, either  physical  or  non-physical.   Examples  of  such  phe- 
nomena include  the  following:   heat  flow  problems  (thermo- 
dynamics), simple  electrical  circuits  (electrical  engineering), 
force  problems  (mechanics) ,  rate  of  bacterial  growth  (bio- 
logical science) ,  rate  of  decomposition  of  radioactive 
material  (atomic  physics) ,  crystallization  rate  of  a 
chemical  compound  (chemistry),  and  rate  of  population  growth 
(statistics) . 

Most  scientists,  engineers  or  applied  mathematicians, 
who  have  faced  the  problem  of  solving  an  ordinary  differential 
equation  numerically,  are  probably  aware  of  the  multitude 
of  techniques  available  for  such  a  problem.  The  abundant 
literature  on  the  subject  of  numerical  solution  of  ordinary 
differential  equations  is  on  the  one  hand,  a  result  of  the 
tremendous  variety  of  actual  systems  in  the  physical  and 
biological  sciences  and  engineering  disciplines  that  are 
described  by  ordinary  differential  equations  and,  on  the 
other  hand,  a  result  of  the  fact  that  the  subject  is  cur- 
rently active. 

The  existence  of  a  large  number  of  methods,  each  having 
special  advantages,  has  been  a  source  of  confusion  as  to 
what  methods  are  best  for  certain  classes  of  problems.   It 
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is  this  observation  which  particularly  motivated  the  writing 
of  this  paper;  thus  it  is  attempted  to  bring  together  the 
well-known  predictor-corrector  methods,  which  form  a  class 
of  numerical  integration  methods  for  solving  ordinary  dif- 
ferential equations,  in  a  consistent  and  comprehensive 
framework. 

This  delimitation  of  the  study  to  predictor-corrector 
methods  is  not  without  basis.   Further  research  in  the  field 
showed  that  the  predictor-corrector  forms  are  the  most 
efficient  among  the  known  integration  methods  in  terms  of 
speed  and  accuracy.   Collectively,  as  a  class  of  integration 
methods,  the  predictor-corrector  sets  are  the  best,  but 
individually  as  a  predictor-corrector  set,  the  choice  for 
the  best  method  varies  depending  on  the  application.   This 
paper  atempts  to  compare  the  overall  efficiency  of  all  the 
well-known  and  most  popular  predictor-corrector  sets  by 
using  a  standard  mode  of  application,  the  same  series  of  step 
size  values,  and  a  set  of  test  ODE '  s  with  many  unusual  and 
interesting  features.   To  introduce  the  numerical  experimen- 
tation of  the  different  predictor-corrector  methods,  a 
comprehensive  analysis  of  each  P-C  set  starting  with  its 
derivation,  error  propagation,  and  stability  is  presented 
in  detail. 

Numerical  experiments  are  conducted  using  the  IBM  360/67 
digital  computer.   The  tremendous  computational  capability 
and  speed  of  this  computer  offered  an  indispensable  tool  in 
conducting  the  experiment  on  a  wide  variety  of  test  ODEs. 
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Single  precision  has  seven  digit  accuracy  and  double  pre- 
cision has  fourteen  digit  accuracy  available  for  computation 

The  paper  starts  with  the  description  of  the  nature  of 
ODE  initial  value  problem.   Then  the  various  P-C  methods  for 
obtaining  numerical  solutions  of  the  problem  are  enumerated. 
To  lead  into  the  derivation  of  formulas  a  brief  review  of 
backward  difference  operator  and  the  well-known  Newton  Back- 
ward Formula  is  presented.   In  turn  the  different  type  of 
P-C  methods  are  discussed  followed  by  the  analysis  of  error 
propagation  and  numerical  stability.   The  stability  bounds 
of  the  P-C  methods  were  established  through  numerical 
experimentation.   The  next  step  is  to  subject  each  P-C 
method  to  a  wide  variety  of  test  ODE.   Then  systematic  com- 
parison of  the  performance  of  each  P-C  method  is  made.   Con- 
clusions are  derived  and  recommendations  are  made  based  on 
the  analysis  of  numerical  results  obtained.   Flow  charts  and 
computer  programs  are  attached. 
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II.   NATURE  OF  THE  PROBLEM 

A.   THE  NUMERICAL  PROBLEM  AND  NOMENCLATURE 

A  linear  first-order  ordinary  differential  equation 
(ODE)  of  the  form 

dy/dx  =  y'(x)  =  f(x,y)  (1-1) 

with  the  initial  condition 

y(x0)  -  y0  (1-2) 

where  f(x,y)  indicates  a  dif ferentiable  function  of  the 
variables  x  and  y  is  commonly  referred  to  as  an  initial 
value  problem. 

A  typical  elementary  differential  equations  text  presents 
several  general  classes  of  methods  for  solving  a  linear 
first-order  ODE.   The  principal  classes  of  methods  are 
(1)  variables-separable,  or  reduction  thereto;  (2)  exact 
equations,  or  reduction  thereto;  and  (3)  solution  by  infinite 
series.   The  student  is  taught  to  apply  the  general  method 
that  appears  best  for  the  solution  of  the  particular  ODE. 
For  example,  the  linear  first-order  differential  equation 

dy/dx  =  xy  (1-3) 

can  easily  be  solved  by  the  variables-separable  method. 
This  is  accomplished  by  rewriting  the  equation  in  the  form 

dy/y  =  xdx 

and  integrating  both  sides  to  obtain 
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lny  =  x  /2  +  c 

where  c  is  an  arbitrary  constant  of  integration, 
solution  can  then  be  written  as 


The  general 


y  =  c.e 


x  /2 


where  c,  =  e  .   The  general  solution  of  such  a  linear  first- 
order  ODE  consists  of  a  family  curves  called  the  integral 

x  /  2 
curves.   The  family  of  integral  curves  y  =  c,e     that  con- 
stitute the  solution  of  equation  (1-3)  is  illustrated  in 
Fig.  1. 


Figure  1.   Solution  Curves  of  y'  =  xy. 


For  each  positive  value  of  c,  a  particular  member  of  this 
family  of  curves  is  determined. 

A  particular  solution  of  equation  (1-3)  can  be  deter- 
mined if  a  condition  on  the  solution  curve  is  specified. 
For  example,  if  it  is  required  that  the  particular  solution 
curve  pass  through  the  point  (0,1),  given  by  the  initial 
condition 


y(0)  =  yn  =  l, 
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x2/2 
then  the  particular  solution  y  =  e     is  obtained  since 

ci  ■  l- 

However,  in  real-life  problems,  many  differential  equa- 
tions that  are  encountered  cannot  be  solved  by  elementary 
classical  methods.   In  such  instances,  one  must  resort  to 
numerical  methods  for  obtaining  one  or  more  particular 
solutions  of  the  initial-value  problems.   A  multitude  of 
techniques  are  available  for  solving  such  a  problem  numeri- 
cally.  This  paper  considers  only  numerical  methods  using 
predictor-corrector  equations.   Even  with  this  restriction, 
a  large  number  of  methods  are  in  existence.   Literature 
research  in  this  subject,  however,  reveals  that  the  most 
popular  and  efficient  methods  are  the  following: 

1.  Euler  Predictor-Corrector  Method 

2.  Milne  Predictor-Corrector  Method 

3.  Nystrom  Midpoint  Predictor-Euler  Corrector  Method 

4.  Hermite  Predictor-Milne  Corrector  Method 

5.  Milne  Predictor-Hamming  Corrector  Method 

6.  Adams  Predictor-Corrector  Methods 

These  methods,  each  having  special  advantages,  have  been  a 
source  of  confusion  as  to  what  methods  are  best  for  certain 
classes  of  problems.   This  paper  will  attempt  to  present  a 
comparative  analysis  of  these  various  predictor-corrector 
methods . 

In  attempting  to  compare  the  various  predictor-corrector 
methods,  such  questions  as  the  importance  of  the  number  of 
function  evaluations,  the  step  size  to  be  used  in  a  numerical 
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calculation,  the  computing  time  required  for  each  method, 
the  influence  of  truncation  and  round-off  error  and  the 
stability  of  various  predictor-corrector  modes  of  computa 
tion  will  be  considered. 
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III.   DERIVATION  OF  PREDICTOR-CORRECTOR  EQUATIONS 

A.   BACKWARD  DIFFERENCE  OPERATOR 

To  lead  into  the  development  of  certain  equations  of 
major  importance  in  numerically  solving  ODE,  a  linear  back- 
ward difference  operator  V  is  defined  by  the  following 
equations 

Vy   =  y   -  y   i  (2-1) 

7n    7n   7n-l  ^    ' 

V2y   =  Vy   -  Vy   ,  =  y   -  2y   ,  +  y  n 
7n    7n    7n-l    7n    7n-l    7n-2 

V3y   =  v2y   -  V2y   ,  =  y   -  3y    +  3y   0  -  y   , 
7n     7n     7n-l    7n     n-1    7n-2    7n-3 


vq    =  vq-l     _  vq-l 

7  r        7n        7n-l 

In  general,  these  formulas  start  with  a  given  sequence  of 

the  y(x  )  =  y   the  ordinates  at  equally  spaced  x   =  x^  +  nh, 
J  v   nJ         J  n  -i     /   r       no 

h  is  a  constant. 


B.   NEWTON  BACKWARD  DIFFERENCE  FORMULA 

Using  these  operators  the  well  known  Newton  backward 
difference  interpolation  formula  is : 

2 

n      a(a+l)V  r~    0-, 

y  =   y   +  aVy   +  — * — r-H y   +  (2-2) 

7n+a   7n     7n      2!      n  K        J 

J  a(a+l) . . .  (a+n- 1)  „n  ,          an 

+  — 5, J- ^ J-   V  y  ,   where  a  =  
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This  formula  is  obtained  by  fitting  a  polynomial  to  the  set 
of  points  (x  ,y  )       at  equally  spaced  {x  }.  (Actually 
the  formula  is  a  polyn  mial  in  a  because  of  the  change  of 
variable.)   Fitting  is  taken  to  mean  that  the  polynomial 
passes  through  the  points  (x  ,y  ) .   Since  there  are 
(n+1)  points  (y  ,y,  ,...,y  ),  this  means  that  if  the  solution 
function  y(x)  is  a  polynomial  of  degree  n  or  less,  the 
formula  above  will  be  exact.   When  y(x)  is  not  a  polynomial, 
or  is  a  polynomial  of  degree  greater  than  n,  the  formula 
will  not  be  exact  (except  at  the  points  themselves).   In 
other  words  an  error  will  be  made  because  a  finite  rather 
than  an  infinite  series  is  used.   In  general  form  this  error 
will  appear  as 

T   =  C  hn+1y[n+1]CO 

n+ 1 
where  C  is  a  function  of  a  and  y      (£)  is  evaluated  at 

some  £  ,  x   <  £  <  x  .   T   is  commonly  referred  to  as  the 

'   o        n    a  } 

truncation  error. 

C.   INTEGRATION  FORMULAS  FOR  ODE 

Integrating  both  sides  or  equation  (1-1)  between  (xn>yn) 
and  (x  +, ,y  +,)  there  results 

yn+1  =  yn  +  /  n+1  f(x,y)dx  (2-3) 

xn 

xn+l 
=  yn  +  f  y' (x)dx 

xn 


=  y   +  h  /   y '  (a)  da 
n     0 
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Attaching  the  index  sequence  of  n . =  0,1,2,...  to  this  equa- 
tion and  noting  that  y   is  known  at  x  ,  (1-2),  it  becomes 
apparent  that  (2-3)  can  be  used  recursively  to  generate 
y, ty 2»y x* ' ' '    as  l°ng  as  there  is  some  way  to  evaluate  the 
integral.   In  other  words  the  solution  of  the  ODE  is  con- 
verted to  evaluating  an  integral. 

Derivations  of  many  of  the  equations  used  to  represent 
(2-3)  now  follows: 

a)   Finite-difference  table  of  backward  differences. 
Assuming  that  (x1,y1),  (x2 ,y2) . . .  (xR ,yn)  are  given,  and 
computing  y'  =  f(x.,y.)  for  i  =  0,1,... n  the  following  table 
is  constructed. 


x      y ' 

o     J  o 


vy-i 


*i       y{  *2y'2 


Vyl 


X2     y2 


Vqy' 
7n 


x   ,    y'  V2y' 

n-1   7n-l  7n 

Vyn 

x      y ' 
n     7n 

These  values  will  now  be  used  to  continue  the  solution  in 
computing  y   , .   In  general,  by  retaining  q    differences 
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in  the  Newton  Backward  Formula  (NBF)  (2-2),  the  following 
finite  series  results. 

yt       =  y  I   +  aVy.   +  0t^|1)   y2y.   +   ... 

7  n+a   7n     7n     2 !      n 

a(a+l)..a(a+q-l)  yq  ,  f2_4. 

q!  7n  ^  *J 

Using  this  interpolating  polynomial  to  approximate  the  inte- 
grand in  (2-3)  [replace  y'(x)  by  y1   ]  yields:   (complete 
details  of  integration  are  presented  in  McCalla  [Ref.  1]) 

q     i 
y  ._!  =  y   +  h   Z   a.V  y  ,   with  a   =  1 
7n+l    7n      .  ~   l   7n'        o 
i  =  0 

where 

r      a  (a+1)  .  .  .  (a+i- 1)  ,     r    .  .  n        /-->  ^^ 
a.  =  /   — » « a 1-   da,   for  i  >  0       (2-5) 

1    0         i: 

Equation  (2-4)  is  in  extrapolation  form  as  can  be  seen  from 

the  table  of  differences,  where  the  end  point   x  ...  is 

r       n+1 

excluded  from  the  interpolating  points. 

The  error  term  associated  with  truncating  after  the  q 
V  is 

T   =  h^2  f1   a(tt  +  l)...(a  +  q)   [q+2]  (  }  d 


0 


But  since  the  coefficient  of  y  "     does  not  change  in  sign 
in  [0,1  ,  it  is  possible  to  write 

TQ  -  aq+1  hi*2/1**23  ft)  (2-6) 

By  actually  calculating  the  a.  in  (2-5)  one  arrives  at 
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Note  that  by  truncating  this  series  at  the  first  difference 
then  (for  q  =  0) , 

y  xi  =  y    +  hy'    with  t    =  \  h2  y[2]  (t)  ■  (2-7) 

7n+l        7n  7n  a        2  7         vw  v         ' 

which  is  the  well  known  Euler's  formula.   Adopting  the  con- 
vention that 

T   =  T(x,h) 

since  a  is  a  function  of  x  and  h.   Also 

T(x,h)  =  \   h2  y[2]a) 

can  be  written  as 

T(x,h)  =  0(h2) 

2 
indicating  that  the  truncation  error  is  of  order  h  . 

Actually  (2-3)  can  be  generalized  by  rewriting  it  in  the 

form 

1 

y  .Li  =  y     +  h  /   y'    da  (2-8) 

'n+1   'n-r        7n+a  k 

-  r 

where  r  is  any  positive  integer.   For  example,  the  case  r=0 

is  the  one  already  discussed.   Following  the  same  procedure 

as  previously  presented,  the  results  for  r  =  1,3,5  are 

obtained  by  substituting  (2-4)  into  (2-8)  where 


_  r   a(a+l)...(a+i-l)j    ,-    .    no 
•  =  /   — * — r-j-*1 -  da,  for  l  =  1,2, ...q. 


Actually  calculating  a-  for  different  values  of  r  there 
results 
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r  ■  Xi   yn+l  =  Xn-1  +  M2+0V+  I  V2+  I  V3+  H  ^  +      (2"9) 


3      3      90 

]y' 


n 


r  ■  3;   yn+1  -  yn_3  +  h[4-4V+  |  V2+0V3+  14  V4  +       (2-10) 

]y' 


r  =  5;   yn+l  =  yn-5  +  ht6-12V+15V2-9V3+  ||  V4  +  OV5  +  (2-11) 

]y' 

The  error  term  associated  with  truncation  of  these  formulas 
is  more  complicated  to  calculate  because  the  integrand 
changes  sign  in  [-r,l],   For  special  cases  however,  as  in 
the  above  for  r  =  1,3,5,  the  coefficient  of  the  r   dif- 
ference is  always  zero.   Thus  it  is  only  necessary  to  use 
the  r-1  differences  of  y'  in  order  to  calculate  a   result 
whose  truncation  error  is  of  the  order  of  r+2  in  h.   For 
r  =  q  =  1 

y   Al    =    y      i    +    2hy*  (2-12) 

7n+l7n-l  7n  v 

with 

3 
T(x,h)    =    !L  y[3]U)     . 

Equation  (2-12)  is  known  as  the  Nystrom  Midpoint  formula. 
Similarly  for: 

r=3;  q=3;  y  x,  =  y   ,  +  4h[y'-Vy'  +  f  V2y']  (2-13) 

n   '  yn+l    7n-3      7n  Jn        3   7n  v 

with 
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T(x,h)  =  11  h5  y[S]U) 

r  -  5;  q  =  5; 

y    =y   +  6hrv'  -  2Vv'  +  —  V  v'  -—  V3v'  +  —  V4v'l 
7n+l   yn-5     lxn    vyn   2   yn   2   yn   20   ynJ 


with 


T(x,h)  =  ^  h7y[7](C) 


Generally,  formulas  derived  using  the  extrapolation  form  of 

the  NBF  are  referred  to  as  open,  explicit  or  predictor 

equation  because  y  jn  occurs  only  on  the  left  hand  side  of 
n  7n+l  } 

the  equation.   In  other  words,  y  ^,    can  be  calculated 

n  '  ' n+ 1 

directly  from  the  right-hand  side  values. 

b)   Extending  the  table  of  backward  differences  by  one 
line  to  include  x   ,  as  an  interpolating  point  yields 


X 

o 

K 

Vy' 

xl 

y{ 

• 

V2y^ 

• 

• 

• 

• 

• 
• 

vV 

7n 

• 
• 
• 

1 

• 

Vyn 

vV 

7n+ 

X 

n 

y' 

7n 

• 

Vyn+ 

1 

'2y  ♦ 

J  n+ 

1 

n+ 

1 

y' 

/n+ 

1 

NBF  now  in  interpolation  form  yields  a  finite  series  as 
follows : 
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yn+a   yn+l    LOt  iJ  yn+l      2!      yn+l 

(a-1) (a) (a+1) . .. (a+q-2)   q   ,         r?-141 
q!  yn+l       ^  14J 

Repeating  the  procedure  used  in  the  extrapolation  form 
(2-4),  analagous  results  for  r  =  0,1,3  and  5  are 

r  =  0;    y   ,  =  y   +  h[l  -  i  V-  —  V2  -  —  V3-  -i5_  v4  + 
rn+l    yn   "LJ-    2     12      24      720 


r  -  1;   y  Al  =  y  n  +  h[2-2V+  4-  V  -  0V"   ^  v 
'    7n+l    7n-l 


|vW--L. 


...]yn+1  (2-16) 

r  =  3'    yn+l  "  yn-3  +  h[4"8V  +  f  V2  -  §  V3  +  ^  V4  -  OV5 


•••]yA+i  ^2-17^ 

r  =  5;   y  ^,  =  y   .  +  h[6  -  18V  +  27V2  -  24V3  +  i^  V4 
'   7n+l   7n-5     l  10 

Some  interesting  and  useful  results  may  be  obtained  by  trun- 
cating after  a  certain  number  of  differences.   In  the  r  =  0 
case  truncation  after  the  first  difference  (q=l)  yields 

y  a.1  =  y  +  t  [yVi  +  y']  (2-19) 

7n+l    7n    2  l/n+l    7nJ  v     ' 


with 


3 
T(x,h)  =  ^2  y[3]C^) 
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which  is  called  the  modified  Euler  formula.   For  r  =  1,  . 
truncating  after  the  second  difference  (note  that  the  third 
difference  is  zero)  yields 


yn+l   yn-l  +  3 


'n+1    7n   7n-l 


(2-20) 


with 


T(x,h)  = 


—  y[5](C) 

90  y         ^J 


which  is  popularly  known  as  Milne's  corrector.   Formulas 
resulting  from  the  use  of  NBF  interpolation  form  are  referred 
to  as  closed,  implicit  or  corrector  equations  since  y_  +  -i 
occurs  on  both  sides  of  the  equation.   In  other  words 

unknown  y  ,,  cannot  be  calculated  directly  since  it  is  con- 

;  n+1  ' 

tained  within  y'  ... 
7n+l 

c)   In  sections  a  and  b  the  equations  were  obtained 
directly  from  manipulating  interpolation  or  extrapolation 
polynomials  using  NBF.   However,  another  method,  called  the 
method  of  undetermined  coefficients,  applied  to  the  gen- 
eralized linear,  K-step  differential -difference  equation 
with  constant  coefficients 


yn+l 


K 


K 

Z 


a-y    +  h  z   3-y'  . 


(2-21) 


will  yield  all  the  formulas  derived  so  far  and  at  the  same 
time  can  be  used  to  obtain  all  other  predictor-corrector 
equations  by  suitable  choice  of  the  parameters  a-  and  3.. 
But  since  the  previous  formulas  derived  are  sufficient  for 
applications  and  only  the  Hermite  predictor  is  of  interest 
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using  the  present  method,  the  tedious  and  long  process  of 
derivations,  which  can  be  found  in  most  numerical  analysis 
books  [Refs.  1,2,3],  are  left  out  and  only  the  results  for 
the  Hermite  predictor  from  [Ref.  4]  are  presented.   There 
it  was  found  that 

a   =  -4    a   =  5    a   =  a   =  o 
o  1         2     3 

30  -  4    Bl  =  2    6_!  =  B2  =  ---  =  3g  =  0 

yielding 

y    =  -4y   +  5y    +  h  [ 4y '  +  2y '  ,]  (2-22) 

7n+l     7n    7n-l      l/n    7n-l  v     ' 

with 

4 
T(x,h)  =  \  y[4]U) 

d)   It  is  of  interest  to  note  that  the  numerical  methods 
developed  for  solving  the  initial  value  problem  differ  in 
their  requirements  of  the  number  of  starting  values.   For 
example,  Euler's  formula  (2-7)  needs  only  the  initial  con- 
dition y(x  )  =  y   and  h  to  start  the  solution.   Methods  that 
J  v  o    J  o 

determine  y  +, ,  when  only  one  point  (x  ,y  )  and  step  size 
h  are  known,  are  commonly  called  one-step  method. 

On  the  other  hand,  the  Nystrom  Midpoint  formula  (2-12) 
needs  starting  values  (x  _i>yn_i)  an^  (xn>yn)  t0  continue 
the  solution.   Methods  that  require  step  size  h  and  more 
than  one  point  (x  ,y  ) ,  (x  _, ,y  _,)....  in  order  to  compute 

y  ^n  are  called  multi-step  methods. 

'  n+1  r 

It  is  clear  then  that  a  multi-step  method  requires  a 
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one-step  method  to  provide  the  necessary  starting  values 
in  its  computation. 

Literature  research  showed  that  the  Runge-Kutta  method 
is  the  best  among  the  available  one-step  methods.   Hence  the 
fourth  order  Runge-Kutta  method  based  on  the  Taylor's  series 
expansion  truncated  after  terms  of  the  fourth  order  (deri- 
vation of  this  method  is  given  in  Ref.  3)  was  used  as  a 
starting  method  for  the  multi-step  methods  considered.   This 
has  the  form 


>Vl   =   >^n   + 


K1+2K2+2K3+K4 


where 


K-,  =  hf(x  ,y  ) 
1      ^  n,7n^ 


K2  =  hf(xn  *\.  yn  +ji) 
K3  =  hf(xn  ♦  \,  yn  ♦  £) 
K4  =  hf(xn  ♦  h,  yn  ♦  K3) 


to  solve  the  initial  value  problem 


y'  =  f(x,y)   with  y(xQ)  =  yQ 


yielding  as  many  starting  values  as  needed 
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IV.   PREDICTOR- CORRECTOR  METHODS 

A.   DEFINITION 

Algorithms  that  use  both  an  explicit  formula  and  an 

implicit  formula  are  called  predictor-corrector  methods. 

The  solution  approximation  computed  by  the  explicit  formula 

is  denoted  y  , ,  and  is  called  a  predictor.   This  predictor 
7n+l  r  r 

is  then  used  initially  in  the  right  side  of  an  implicit 

formula  that  computes  a  corrector  y   - .   The  implicit 

formula  can  be  repeatedly  applied,  using  y  +.  from  the 

preceding  iteration  on  the  right  side  and  computing  a  new 

y0^,  on  the  left. 
'n+1 

To  illustrate  the  procedure,  the  Euler  P-C  algorithm 
(one-step  method  predictor)  and  the  Nystrom  P-C  algorithm 
(multi-step  method)  are  used. 

Given  the  ODE 

y'  =  f(x,y),  with  initial  solution  point  (x  ,y  ) 

1)   Euler  P-C  method  uses  the  equation  (2-7)  as  predic- 
tor, and  (2-19)  as  corrector 

yn+i  ■  yn  +  I  iyi+i  +  w  (2-24) 

Since  (2-23)  requires  only  one  solution  point  (x  ,y  )  and 
step  size  h  (assumed  to  be  constant  for  specific  application) 
which  can  be  chosen  arbitrarily  (the  choice  actually  affects 
the  stability  and  error  propagation  of  the  method  as  will 
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be  illustrated  later),  y'  =  £(x  ,y  )  can  be  evaluated.   Hence 
(2-23)  can  be  computed  yielding  the  result  y  +1   (p  means 
predicted  value) .   Then  the  derivative  of  y  +1  is  computed 
yielding  the  result  y '  - : 

yn+l  =  f(xn+l'yn+l^  * 
Next  (2-24)  is  used  to  compute  a  corrected  value  called 
y   +,   (c  means  corrected  value) 


?n+\   '   yn  +  7  [K*l   +  yAl  • 


At  this  point  there  are  two  alternatives.   First,  the  compu- 
tation  can  be  terminated  using  y  +,  as  the  true  approxima- 
tion value  for  y  , , .   Then  the  computation  continues  to 

7n+l  r 

determine  the  next  solution  point  y  +?   by  repeating  the 

entire  procedure.   Second,  the  corrected  value,  y  ^,  of 
r  '  ' n+1 

y_  +  -j  may  not  be  acceptable  in  the  sense  that  the  equality 
of  both  sides  of  (2-24)  is  not  satisfied  to  as  many  digits 
as  may  be  desired.   In  other  words,  the  corrector  has  not 
converged  to  the  desired  accuracy  (convergence  of  the 
corrector  formula  is  an  important  topic  and  discussion  of 
this  aspect  is  deferred  to  a  later  section) .   The  next  step 
then  is  to  improve  the  value  y  +,  on  the  left  side  of  (2-24) 
by  taking  the  derivative  of  yC+-i  ,  calling  the  result  y'  ,, 
and  then  calculating  an  improved  value  of  y   , .   This  is 
repeated  until  convergence  occurs  to  as  many  digits  as 
desired.   This  convergence  term,  designated  as  EPS,  is 
tested  against  the  absolute  difference  between  the  last  two 
approximations . 
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2)   The  Nystrom  P-C  method  has  (2-12)  as  predictor  and 
(2-19)  as  corrector. 

>vi  '  yn  +  I  Oa+i  +  *i]  (2-26^ 

It  is  clear  that  before  (2-25)  can  be  evaluated  another  solu- 
tion point  (x   -,  ,y   -, )  is  needed  in  addition  to  the  given 
r      v  n- 1 '' n- 1'  & 

solution  point  (x  ,y  ,)    of  the  ODE.   Another  one  step  method 
r      v  n,7n-l^  r 

is  therefore  needed  to  provide  the  additional  starting 
values.   The  Runge-Kutta  method  can  be  used  to  accomplish 
this.   Once  these  starting  values  are  known,  the  same  pro- 
cedure as  in  Euler  P-C  method  is  followed  to  compute  the 
next  point  using  equation  (2-25)  as  predictor  and  (2-26)  as 
corrector . 

B.  A  SIMPLE  PREDICTOR-CORRECTOR  SET 

The  predictor-corrector  methods  discussed  so  far  do  not 
take  into  account  the  truncation  error  incurred  in  retaining 
r   differences  of  the  infinite  series  of  the  NBF.   The 
application  is  straightforward,  predict  then  correct,  and 
as  such  it  can  be  called  a  simple  predictor-corrector  set. 

C.  MODIFIED  PREDICTOR-CORRECTOR  SET 

The  P-C  set  given  by  (2-12)  and  (2-19)  has  a  local  trun- 
cation error  of  the  same  order  as  the  set  (2-23),  (2-24) 
(T(x,h)  =  0(h3)).   Hamming  [Ref.  4]  has  used  this  feature  to 
expand  the  method  of  calculation.   The  local  truncation 
errors  for  the  predictor  (2-12)  and  corrector  (2-19)  are 
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where 


and 


T(x,h)p  =  Tn+1)p  =  §-  Y[3]  Up)        m  (2-27) 


n-1    p    n+ 


T(x-h'c  =  Vl,c=-^l31Kc)  C2"28^ 


where 


n   ^c    n+1 

It  follows  that  the  exact  value  y(x  .,  )  ,  of  y  at  x  . .  .  is 

'  *•  n+1' '     }  n+1 ' 

given  by 

y(xn+1)  -  yP+1  +^y!3)(cp) 

3 
/-     -\     c      h    [  3]  • ,.  n 

y(xn+1)  ■  yn+1  -  -yyy1  (cc) 

From  these  equations,  one  obtains 

c  h3      [3]  rr    .  p  h3      [3]  rr    , 

yn+l    "    12   y         (y    =   yn+l    +    3~y         (S} 


or 


&i  -Crry[J1(S'  +  T^y[3!^         t2"29) 

Now  assuming  that 

y[31(Cp)  =  y[31(Cc)  =yI31(C) 

Over  the  interval  of  interest  (i.e.,  the  third  derivative 
does  not  vary  greatly  over  the  interval),  (2-29)  becomes 


31 


where 

n-1        n+1 

Even  at  this  point  in  the  development  some  interesting  fea- 
tures can  be  identified.   First  one  notices  that  (2-30)  can 
be  obtained  only  in  case  the  order  of  the  predictor  equation 
is  equal  to  the  order  of  the  corrector  equation.   Second, 
there  is  now  a  measure  of  the  local  truncation  error  in  terms 

of  y  .,  and  y^  ,  which  are  explicitly  calculated  in  the  P-C 
'n+1     J n+1  r  J 

algorithm.   This  may  be  seen  more  explicitly  by  comparing 
(2-27)  with  (2-30). 

"  f  (f  h3y  3  (£)) 

O    rp 

4  n+l,p 
Similarly  comparing  (2-28)  with  (2-30) 

Since  the  local  truncation  error  can  now  be  estimated,  the 
next  question  is  how  to  use  this  information  to  improve  the 
predicted  and  corrected  values.   Considering  the  predicted 
value  first,  it  can  be  seen  from  (2-31)  that 

T       =  1  (yc    -  yP   ) 
n+l,p    5  ^7n+l    7n+lJ 

Assuming  that  the  local  truncation  error  remains  approxi- 
mately constant  over  two  steps,  then 
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n+l,p    5  lyn+l  "  yn+lJ    5  l7n  "  ynj 

The  predicted  value  can  then  be  improved  or  modified  by 

adding  the  term  ■=■  (y      -    y*,)  to  yP  n  using  only  information 

b  5   n   J  n  '  n+1     b  J 

calculated  previously. 

The  corrector  can  be  modified  in  essentially  the  same 
manner  since 

T       =  -  I  fyc    -  yp  ) 
n+l,c     5  Lyn+1   yn+lJ 

can  be  added  to  y   -  to  improve  this  value.   On  this  basis 
the  overall  P-C  calculation  can  be  written  in  the  following 
algorithmic  steps: 

Predict:      p  .,  =  y   n  +  2hy' 
*n+l    7n-l     7n 

Modify:       m  ^n  =  p  ^n  -  ■=•  fp  -zj) 
7  n+1   *n+l    5  Vin  n' 

Reevaluate    m1    _  ,,      „    >> 
Derivative:     n+1     v  n+1'   n+1' 


h 
'n+1    7n    2 


Correct:      c  ...  =  y   +  %   (m'  ,  +  fn) 

7n    2    n+1      J 


Modify:       yC.i  =  c  .-  +  —  (p  .-  -  c  _,, ) 
7        7n+l    n+1    5   *n+l    n+1' 

Then  iterate  the  corrector  to  convergence.   To  simplify  the 

presentation,  the  symbols  m  , , ,  p  ,.  and  c  ,,  have  been  used 
r  '      3  n+1'  rn+l      n+1 

The  use  of  the  modified  predictor-corrector  set  seems 
to  be  very  attractive  but  it  must  be  noted  that  this  method 
has  two  limitations.   First,  it  rests  on  the  fact  that  the 
original  P-C  equations  have  equal-order  truncation  errors, 
and  second,  it  neglects  round-off  error. 
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D.   HAMMING  MODIFIED  PREDICTOR-CORRECTOR  SET 

Hamming  used  equations  (2-10)  and  (2-16),  which  are  known 
as  the  Milne  P-C  set  when  both  are  truncated  after  the  third 
differences,  and  he  revised  the  corrector  to  remove  the  severe 
stability  problems,  which  will  be  seen  later. 

To  develop  the  revised  corrector,  Hamming  started  with 
the  generalized  equation  (7-21)  and  obtained 

>Vl    =    %yn   +    al>Vl    +    °2>V2    +   hlfi-lK+l   +    VA 

Using  the  method  of  undetermined  coefficients  [complete  deri- 
vations can  be  found  in  Ref .  4] ,  he  obtained 

%  "  J  (9-9V        B-i  "  h  (9"V 

ai  =  ai  eo  =  n  C9  +  7ai} 

a2    -    -    \   (1-c^)         g1      =   i    (-9    +    17aj) 

with 

(-9   +    Sa.)       5    r5, 

T(x,h)  =      360       hys'm 

with  a,  as  a  free  parameter.   Hamming  then  considered  the 

values  of  a  ,a0,...$,  as  functions  of  a,  for  a,  =  1,  9/17, 
o '  2 '     1  1      1 

1/9,  0,  -1/7,  -9/31,  and  -6/10.   After  looking  at  the  result- 
ing coefficients  in  terms  of  equal  magnitude,  number  of 
zeros,  etc.,  the  case  ou  =  0  was  chosen  as  the  best  value 
since  it  yields  a  greater  region  of  real  negative  stability 
and  a  wide  range  of  relative  stability.   It  can  be  verified 
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that  the  case  a,  =  1  yields  the  Milnes  corrector  (2-10) 
truncated  after  the  third  difference.   With  a-,  =  0  the 
Hamming  predictor-corrector  set  is  of  the  form 

which  is  (2-10)  in  the  case  r  =  3,  q  =  3,  of  the  NBF  extra- 
polation form,  with 

T(x,h)  =  1|  h5y[5](C) 
and 

yCxi    =   J    [9y    -y      ,]    +   |  h    [y»         +    2y'-y!    ,]  (2-34) 

7n+l        8    L    7n   7n-2  8         l/n+l  7n   7n-lJ  v  J 

with 

T(x,h)   =  ^  y[5](C)    . 

Next,  using  the  procedure  discussed  in  the  previous  section 
on  the  modified  predictor-corrector  set,  Hamming  develops 
his  modified  predictor-corrector  algorithm  as  follows: 

Predict:     p  Al-  =  y   _  +  4  h  [2y'-y'  -.  +  2y '  ,]    (2-35) 
*n+l    7n-3    3    l/n/n-l   7n-2     v     * 

112 
Modify:      m  _,  =  p  .-  -  (t-otO  (p  -c  ) 
J  n+1    rn+l    ^1217  K^n      nJ 

Reevaluate     ,      r,  > 

Derivative:   mn+l     lxn+l,mn+lJ 

Correct:      c  ±1  =  ~  [9y  -y   -]  +  -J  h  [m'   +2y'-y'  ,] 
n+1    8    7n  7n-2     8     n+1   7n  7n-l 

c  9 

Modify:      y   ..  =  c   ,  +  r—r    [p   n-c   ,1 
7       7n+l    n+1    121   Fn+1   n+lJ 

Then  the  corrector  may  be  iterated  to  convergence  as  desired 
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E.   P-C  SETS  CONSIDERED  IN  THE  NUMERICAL  EXPERIMENTS 

Having  set  forth  the  necessary  foundation  for  the  pre- 
dictor-corrector equations,  a  list  of  the  P-C  sets,  which 
were  analyzed  and  actually  applied  to  the  numerical  solution 
of  the  initial  value  problem  in  (1-1)  is  now  summarized.   In 

each  case  the  operator  is  evaluated  in  terms  of  y'  -,  ,  y', 

r  ' n+1 '    J n* 

yn-l'*  * ' 

P-C-I:      Euler  P-C   set. 

The   predictor    is    equation    (2-7)    and   the    corrector   is 
(2-24), 

yn+l   =   yn   +   hyn    ;  T(x'h)    =   1  hV21  CO 

P-C-II:   Milne  P-C  set. 

Equation  (2-10),  in  the  case  r  =  3,  q  =  3  of  the  NBF 
extrapolation  form,  is  used  as  the  predictor,  and  the  cor- 
rector is  (2-16)  where  r  =  1,  q  =  3,  of  the  NBF  interpolation 
form.   This  yields  the  equations 

P-C-III:   Nystrom  Predictor-Euler  Corrector  Set. 

The  Predictor  is  (2-12)  and  the  corrector  is  (2-24)  which 
yields 

Ci  -  yn-i +  2hyA  ■    T^^  -  r  y[31(« 
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.P-C-IV:   Hermite  Predictor-Milne  Corrector. 

The  predictor  used  is  (2-22)  and  the  corrector  is  (2-16) 
yielding 

yl*i  =  •%  +  5>Vi  +  h  "/A  +  ^A-i1 

T(x,h)    =  t-  y[4!ca 

0 

yn+l    "   >Vl    +   I    t^A-1    +    4^n   +   W 

T(x.h)   =  ±-  y[5](£) 

P-C-V:   Hamming  Modified  Predictor-Corrector  Set. 

Equation  (2-33)  is  the  predictor  and  (2-34)  is  the  cor- 
rector. 

*S+i  ■  >V3  +  r  [2yA  "  ri-i  +  2^-2^ 

T(x,h)    =i|h5y[51(0 

>Vl    "   1    t^n^n-2^    +   f    ^A+l    +    2^n    "    Y^-l* 

T(x,h)    =   ^.  yt5](0 

Adams  Form  (Adams -Bashforth  Predictor-Adams -Moulton  Corrector 
Set)  . 

The  A-B  formulas  are  (in  the  case  r  =  0)  equations  of 
the  NBF  extrapolation  form  and  the  A-M  formulas  are  (in  the 
case  r=0)  equations  of  the  NBF  interpolation  form.   Three 
methods  are  considered:   the  second,  third,  and  fourth  order 
A-B  predictor  and  A-M  corrector. 
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P-C-VI:   Second  Order  Adams  P-C  Set. 

The  predictor  and  corrector  are  obtained  in  the  case 
q=l  yielding 

T(x,h)  =  JL   h3y[3)(C) 


yn+i  =  yn  +  I  I7i+1  ♦  /A) 

3 

T(x,h)  =  -i-  y[31(0 
12 


P-C-VII:   Third  Order  Adams  P-C  Set. 

The  case  q  =  2  in  both  the  extrapolation  and  interpola- 
tion of  the  NBF  form  are  used  as  the  predictor  and  corrector 
equations  respectively,  resulting  in  the  formulas 

T(x,h]  -  ^h4yI41(0 

4 
T(x,h)  =  ^  y[4]  CO 

P-C-VIII:   Fourth  Order  Adams  P-C  Set. 

The  case  q  =  3  in  both  NBF  forms  yields  the  predictor- 
corrector  equations 

yn+i  ■  yn  +  ta  [55yA  -  59^-i +  "yA-z  -  ^i-s1 
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T(x,h)  -  ||i  h5y  5  (O 


^n+l  "  ^n  +  £  ^A+l  +  19>^n  "  5^-l  +  Xi-23 


T(x,h)  =  li|  h5y  5  CO 


From  the  P-C  sets  above  an  interesting  fact  is  observed. 
Namely,  P-C-I,  P-C-III,  and  P-C-IV  used  different  predictor 
equations  but  have  the  same  corrector  equation  while  P-C-II 
and  P-C-V  used  the  same  predictor  equation  but  have  dif- 
ferent corrector  equations.   The  significance  of  this 
observation  will  be  clearly  demonstrated  later  on. 
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V.   NUMERICAL  STABILITY  OF  PREDICTOR-CORRECTOR  METHODS 

In  the  numerical  solution  of  an  ODE,  a  sequence  of 

approximations  y   to  the  true  solution  y (x  )  is  generated. 

Roughly  speaking,  the  stability  of  a  numerical  method  refers 

to  the  behavior  of  the  difference  or  error,  y   -  y(x  ),  as 

'  7n    7  K   nJ  ' 

n  becomes  large.  In  order  to  begin  the  discussion,  the 
various  types  of  errors  which  are  incurred  in  numerical 
integration  of  (2-4)  must  be  considered. 

The  errors  incurred  in  a  single  integration  step  are 
of  two  types: 

1.  The  local  truncation  error  or  discretization 
error  -  the  error  introduced  by  the  approximation  of  the 
differential  equation  by  a  difference  equation. 

2.  Errors  due  to  the  deviation  of  the  numerical  solu- 
tion from  the  exact  theoretical  solution  of  the  difference 
equation.   Included  in  this  class  are  round-off  errors,  due 
to  the  inability  of  evaluating  real  numbers  with  infinite 
precision  with  the  use  of  computer  (i.e.,  computers  usually 
operate  on  fixed  word  length) ,  and  errors  which  are  incurred 
if  the  difference  equation  is  implicit  and  is  not  solved 
exactly  at  each  step. 

If  a  multi-step  method  is  used,  an  additional  source  of 
error  results  from  the  use  of  an  auxiliary  method  (usually 
a  single-step  method,  e.g.,  Runge-Kutta  method),  to  develop 
the  needed  starting  values  for  the  multi-step  method. 
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A.   PROPAGATION  OF  ERROR  IN  A  ONE-STEP  METHOD 

The  propagation  of  error  in  a  one-step  method  can  be 
analyzed  by  studying  the  particular  representative  one-step 
method  of  Euler  (2-7). 

y  A-  =  y   +  hf(x  ,y  )  (2-36) 

7n+l    7n      v  n,7n7  *•     ■* 

This  can  be  done  by  determining  the  relation  of  the  error 

at  step  n+1  to  the  error  at  step  n.   To  do  this  let  Y 
v  r  n 

denote  the  true  solution  of  the  initial  value  problem 
y'(x)  =  f(x,y)  with  y (x  )  -  y  .   Then  the  total  accumulated 
solution  error  e   at  step  n  is  defined  by 

The  numerical  values  computed  by  Euler's  algorithm  (2-36) 
satisfy  the  relation 

y  Al  =  y   +  hf(x  ,y  )  -  R  Al  (2-38) 

7n+l    7n      v  n,7n7     n+1  v     J 

where  R  ^,    denotes  the  round-off  errors  resulting  from 
n+1  to 

evaluating  (2-36).   Similarly  the  true  solution  values 
satisfy  the  relation 

Y  Al  =  Y   +  hf(x  ,Y  )  +  T  Al  (2-39) 

n+1    n      v  n*  n7     n+1  v     7 

where  T  .-  denotes  the  local  truncation  error  in  (2-36). 
n+1  v     * 

Subtracting  (2-38)  from  (2-39)  yields 

Y„xi-  7  .n  =  Y   -  y  +  h[£(x  ,Y  )-  f(x  ,y  )]+E  ... 
n+1   7n+l    n   7n     L  *"  n'  nJ         K   n'7n7    n+1 

(2-40) 

where  E  .-  =  T  A-  +  R  .,, .   Inserting  (2-37)  in  (2-40)  yields 

n+1    n+1    n+1  6  v     7     v     J    J 

E  x1  =  £   +  h  [f(x  ,Y  )  -  f(x  ,y  )]  +  E  x1     (2-41) 
n+1    n      l  v  n'  n7     *■  n'7n7     n+1     v 
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Applying  the  mean  value  theorem  to  the  terms  inside  the 
bracket  of  (2-41),  this  relation  between  successive  errors 
can  be  written  as 

E  A,  =  I      +  h  [  (Y  -y  )  £  (x  ,y  )]  +  E  x-       (2-42) 
n+1    n     L  *■  n  ' nJ      yK  n'7nJJ  n+1      v    * 

where  f   denotes  3f/8y  and  y   lies  between  y   and  Y  .   If 
y  '    ;  n  7n      n 

|f  (x,y) I  <  C  and  I  E  ^,  I  <  E  where  C  and  E  are  both  positive 
1  y v  >J J  '  —       '  n+1 '  —  y 

constants,  the  error  expression  above  can  be  replaced  by  a 
related  first-order  difference  equation  (Henrici   Ref.  5 
gives  an  excellent  treatment  of  the  solution  of  difference 
equations) . 

e  x1  =  e   [1+hC]  +  E  (2-43) 

n+1    n  K  J 

Under  the  stated  assumptions  it  follows  from  (2-42)  that 

|Z  _,,  I  <  |S  I  [1  +  hC]  +  E 
i  n+i  i  __  i  ni  i     j  . 

Now  if  I Z  I  <  e   it  follows  that 
1  n '  —  n 

IW  £  en  [1  +  hC]  +  E 

Therefore,  from  (2-43) 

|E  •  , I  <  e  .-  (2-44) 

1  n+1 '  —  n+1  v     J 

It  follows  then  from  (2-44)  that 

|£  I  <  e  ■*■    1 1,1  <  en  ■*■     ....where  the  symbol  ■*■   means 
i  o '  —   o    '  1 '  —   1         -T-    -ui 

implies  that. 

That  is,  the  propagated  error  is  bounded  by  the  solution  of 

the  related  first-order  difference  equation.   Now  since 

Y  -y^  =  £   =  0,  the  condition  that  |Z  I  <  e   will  be  satis- 
o  '  o    o     '  '  o '  —  o 

fied  by  setting  e   =0,  which  is  the  initial  condition  for 

the  difference  equation. 
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Defining  G  =  1+hC,  the  difference  equation  (2-43)  can 
be  written  in  the  form 

e  ,,  =  Ge   +  E   with  initial  condition  e   =0 
n+1     n  o 

The  solution  of  this  first-order  difference  equation  can 

be  found  by  successive  substitution  as  follows 

e,  =  Ge   +  E  =  E 
1     o 

e2  =  Ge1  +  E  =  (G+1)E 
e3  =  Ge2  +  E  =  (G2+G+1)E 

e   =  Ge   ,  +  E  =  (Gn_1+  Gn_2  +  ...+  G+1)E 
n     n-1       v  ' 

The  solution  e   can  be  seen  to  be  a  geometric  progression 
and  as  such  can  be  written  as 


(2-45) 

It  follows  then  that  the  propagated  error  in  Euler's  algor- 
ithm (2-4)  is  bounded  by  the  expression 

|£    |  <  f(l+hC)n+1-l 
|Ln+ll  -  I     hC 

We  shall  illustrate  this  estimate  by  a  numerical  example. 
Given  the  initial  value  problem 

y'  =  f(x,y)  =  -y  with   y(0)  =  1 

we  apply  Euler's  formula  (2-36).   It  is  required  to  express 

the  error  at  x.,  =  1  as  a  function  of  h  assuming  that  the 

-  7 
round-off  error  |R  +, |  £  10    for  single  precision  (double 

precision  yields  IR  x1|  <  10~   )  in  the  IBM  360/67  computer 

r         /       i  n+1 I  —      J  r 
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Computing 

C  =  max  | £  (x,y)| 

yields 

C  =  1 

The  truncation  error  associated  with  (2-36)  is  given  by 

T  xi  =  "  t  h2y"  (€  )   where  x  <  £   <  x  J, 
n+1     2   7   ^n'  n  —  ^n  —  n+1 

Knowledge  of  the  second  derivative  y"(x)  is  needed  to  eval 
uate  this.   Thus 

y"(x)  =  a|  (y'(x))  =  £   (f(x,y)) 

=  fx  +  £y  dy  . 
dx 

Since  f   =  0  and  f   =  -1,  this  reduces  to 
x  y      ' 

y"(x)  =  y. 
Using  (2-45) ,  one  obtains 


<  (1+h)  -1   in-7    1  ,2      I 

=  - f 10    +  o-  h  max   y 

n       h  2        ' J 


0<x<l 


To  evaluate  max | y | ,  it  is  seen  that  the  analytical  solution 
of  the  initial  value  problem  is 

-x 

y  =  e 

Then  for  x  =  1 

y  -  e"1  -I 

'  e 
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It  is  clear  that  max  |y|  is  attained  at  x  =  0 ;  that  is 

max  |y|  =  1. 
Substituting  into  the  original  expression  yields 


,  <  [10-7  +  lv](Ili2li^l, 

n    1         2  ]\  h    j 


(2-46) 


Rewriting  this  equation  yields 


•nM^*!!"!*)11-!) 


(2-47) 


Thus  the  propagated  error  incurred  in  stepping  forward  from 
x  =  0  to  x  =  1,  as  a  function  of  h,  is  bounded  by  the  above 
expression. 

The  revised  form  of  the  propagated  error  bound  shows  the 
influence  the  various  errors  have  in  the  numerical  solution 
of  the  ODE. 

Rewriting  (2-47)  in  the  form 


R 


n 


T1  +  TA+i)^1+h)n- 


(2-48) 


where   R  , ,  is  still  the  round-off  error 
n+1 

T'+,  is  the  truncation  error  reduced  by  0(h) 

and  ignoring  the  contents  of  the  second  parenthesis,  the 
equation  suggests  that  as  h  -*■   0 ,  the  truncation  error 
decreases  toward  zero  whereas  the  round-off  error  increases 
toward  infinity.   To  minimize  the  error,  it  is  clear  that  h 
must  be  chosen  so  that  the  truncation  error  and  the  round- 
off error  have  equal  orders  of  magnitude.   A  numerical 
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experiment  was  conducted  to  verify  this  idea.   Finally  it 
can  be  seen  that  the  accumulated  error  is  not  simply  equal 
to  the  sum  of  the  local  truncation  error  and  round-off  error 
but  must  be  computed  as  given  by  (2-48)  which  is  the  solu- 
tion of  the  first  order  difference  equation  corresponding  to 
(2-36) . 

B.   PROPAGATION  OF  ERROR  IN  THE  MULTI-STEP  METHOD 

As  shown  in  equation  (2-21)  the  most  general  representa- 
tion of  the  Predictor-corrector  methods  are  of  the  form 

K  K 

y    .-    =      E      a.y       .    +   h      E        B-y'     .  (2-49) 

7n+l         .    n      n/n-i  .      ,      i7n-i  *■  J 

i=0  i=-l 

However  by  confining  the  study  to  a  form  represented  by 

K 

y  x1  =  y    +  h   Z    B-f   .  ,   f   .  =  y'  .     (2-50) 
7n+l    7n-r      ._  ,   l  n- i  '    n-i   7n-i     v     J 

simplicity  is  obtained  without  loss  of  generality.   It  can 
be  verified  that  most  of  the  formulas  in  Section  III  are  in 
the  form  (2-50).   For  example,  the  case  r  =  0  and  3  _ -.  =  0 
the  Adams -Bashforth  formulas  are  obtained  as 

K 

yJ.1=y   +h   £   B-f 
7n+l    J n      .  0  Mi  n- i 

while  in  the  case  r  =  0  and  3_-i  ¥    0  the  Adams-Moulton  forms 
result. 

The  propagation  of  error  in  the  multi-step  method  can 
then  be  analyzed  by  studying  equation  (2-50).   Following  the 
same  definition  as  in  the  previous  section  with  regard  to 
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Y  ,  y  ,  E  ,  R^,,  T^,,  E  ^,  ,  C,  and  E, 
n*  /n'   n'   n+1*   n+1'   n+1'   '       ' 

repetition  can  be  avoided. 

The  solution  values  generated  by  (2-50)  satisfy  the 

relation 

K 
>Vl  =  >Vr  +  ^-l^n+l'W  +  h  ±lQ   M^n-i^n-i.) 

-  Rn+r  (2"51) 

Similarly  the  solution  values  Y   satisfy  the  relation 

'  n  J 

K 
Vl    =   Vr   +    hS-l£K+l>W  +  h    ±lQ    MK-i'Vi7 

+  Tn+r  (2-52) 

Subtracting  (2-51)  from  (2-52)  yields  the  equation 

Y  Ai  "  y  Al  ■  Y    -  y    +  h$  -,(f(x  .-  ,Y  .- ) 
n+1   7n+l    n-r  7n-r     -1^  v  n+1'  n+17 

K 

"    f (x   x1  ,y   ,,))    +    h      E      6- (f (x      • ,Y      .) 
v   n+l,7n+l77  .    n      iv    v   n-i'    n-i7 

i  =  0 

-  f(x   •)  +  T  x1  +  R  .-  (2-53) 

*•  n-i7     n+1    n+1         v     J 

By  application  of  the  mean  value  theorem  to  (2-53)  and  using 

the  definition  of  E   and  E  J n  one  obtains 

n  n+1 

£    .-    =    E  +   h$    -  (Y   ^,-y   ^-,)fy(x   A-  ,y   ±1)  (2-54) 

n+1  n-r  H-lv   n+1    7n+l7    3K   n+l,7n+l7  ^  7 

K 

+   h  E      3-  (Y      .  -y      -)fy(x      •  ,y      •)  +  E   .,-, 
.    q      iv   n-i   7n-i7    ,y   n-i,7n-i7         n+1 

Using  the  definition  of  C,E,  and  E  ,  (2-54)  reduces  to 

K 

E  ..  =  E     +  h3  nE  ±1C  +  h   E   g.    .C  +  E 
n+1    n-r     -1  n+1       ._n   i  n-i 
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or 

K 

E  ,-t-  h3  i^  ,iC  =  I        +  h   Z   3-Z   .C  +  E      (2-55) 
n+1     -1  n+1     n-r     .  Q   i  n-i  v     ' 

The  related  difference  equation  to  (2-55)  can  be  written  as 

K 

e  ^(l-hlg  ,|C)  =  e    +  hC  Z  Ig.le   -+E     (2-56) 
n+1^    '  —  1 '  '    n-r      i  =  o     n_1 

If  |Z   . I  <  e   .  (i  =  0,1,... K)  where  K  >  r  it  follows  that 
1  n-i '  —  n-i  v      '  '     '  — 

K 

|Z    A,  I  (    1-hB    ,C|)    <    |Z         |+hC      Z       I  3-  I  |Z       .  I  +  I  El  . 
i    n+ii\-  w -i    ij    _   i    n-r1  •    n    '    i'  '    n-i1       '     ' 


i  =  0 


Thus 


K 

|Z   J(|l-hB  ,C|)  <  e    +  hC   Z   |&.|e   .  +  Ie|. 
i  n+iivi    w -i  \j   _  n-r  i^i  n-!    i  i 

From  (2-56)  it  follows  that 

|sn+1l  C|i-h3.1c|D  <  en+1  (l-hCla.J)  . 

If  hC I  3 _ x I  <  1,  then 

|2  a.-,  I  <  e  .-  (2-57) 

i  n+11  —  n+1  v     ■* 

That  is,  if  the  magnitude  of  the  propagated  error  is  dom- 
inated by  the  solution  of  the  related  difference  equation  for 
K+l  successive  steps,  namely  |Z.|  <_   e-  (i  =  0,1, ...K)  then 
by  induction,  the  same  is  true  for  all  successive  integer 
values  of  i.   Therefore,  a  bound  for  the  propagated  error  in 
the  multi-step  method  can  be  determined  by  obtaining  a 
solution  of  the  related  linear,  inhomogeneous  difference 
equation  of  (2.55). 
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As  with  ODE,  (2-56)  has  homogeneous  and  particular 
solutions  such  that 

e   =  e  u  +  e 
n    nH    np 

The   homogeneous    solution   is    given   by 

enH  =    (l-hC|3_1|)uK+1-UK"r-    hC(|30|uK 

+    ISjJu^1   +    ...+    |6K|)  (2-58) 

The  particular  solution  can  be  obtained  by  assuming  that 

e   =  -X,    then  substituting  in  (2-56) 

n  °  v     J 

(l-h|3_1|C)(-X)-(-A)-hC(-A) [|30|+|31|+...+|3K|]  -E 

K 
A  [hCl  3-,1+hC  I       I  B-  I  ]  =  E 
1       i=0    x 

K 
XhC   Z    I 3- I  =  E 
i-1    X 


X  = 


K 
hC   S   I  3 


i-1   X 

The  general  solution  of  the  difference  equation  of  (2-56) 
can  be  written  as 

en  =  dlul  +  d2u2+  •■•  +  dK+lUK+l+  IT (2"59) 

hC   E   |  3-  | 
i=-l  1 

where  d, ,d- , . . . d„+,  are  determined  by  the  initial  starting 
values  and  where  u. (i=l , 2 , . . . K+l)  are  the  roots  of  the 
characteristic  equation 

ri-hC|3_1|]uK+1-uK"r-hC(|30|uK+|31|uK"1+...+  |eKl)  =  0 
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In  what  follows  the  assumption  is  made  that  the  roots  are 
real  and  distinct  (complex  roots  and  multiple  roots  are 
discussed  in  detail  in  Ref.  4). 

It  may  be  observed  that  the  homogeneous  part  of  the 
characteristic  equation  for  the  accumulated  error  (2-59)  is 
exactly  the  same  as  in  the  characteristic  equation  for  the 
original  difference  equation  of  (2-50).   The  significance 
of  this  point  will  be  used  in  the  analysis  of  numerical 
stability  following  a  numerical  example. 

To  illustrate  the  previous  analysis  consider  the  case 
where  r  =  0  and  K  =  0  in  (2-50)  which  yields 


yn+l  =  >^n  +  2  [£n+l  +  fn] 


(2-60) 


This  is  the  modified  Euler  formula  wherein  y  is  given  by 
the  initial  condition  and  the  needed  additional  point  has 
been  supplied  by  a  predictor  formula. 

Assuming  that  f  (x,y)  =  -C  and  that  E  +,  =  E,  the  related 
difference  equation  yields 

(1  -  |  (-C))  u±   -  1  -  h(-C)  \   =  0 

which  has  root 


ul  = 


(1  -  ffi 


(2-61) 


Using  (2-59)  the  difference  equation  solution  is  then 


e   =  d, 
n    1 


(1  -  ^)n 


(1  + 


hC 
2  -J 


hC  (i  +  i) 
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Analyzing  the  first  term  with  d,  as  constant,  it  remains  to 
show  that  if 

2 

then  the  total  error  e   decreases  with  n  and  convergence  is 

n  ° 

assured.   Furthermore  if  this  condition  is  met  then 

i  -  !£\ 

— - £j<l  (2-62) 

and  the  solution  of  the  difference  equation  decreases  with 
increasing  n.   But  the  expression  on  the  left  of  (2-62)  is 
exactly  the  root  of  the  characteristic  equation  of  (2-60); 
thus  it  is  clear  that  if  u,  <  1,  then  the  solution  of  the 
difference  equation  decreases  with  n.   This  idea  can  be 
extended  to  the  general  characteristic  equation  (2-50)  whose 
solution  may  be  written  in  the  form 

yn  =  dlul  +  d2u2  +  •••  +  dK+luK+l  (2"63) 

It  has  been  shown  earlier  that  the  total  solution  error 

|£  xi  I  <  e  ...i  is  bounded  by  the  solution  of  the  related 
1  n+1 '  —  n+1  ' 

difference  equation.   In  turn  it  was  observed  that  e   depends 
n  n   r 

on  the  solution  of  the  characteristic  equation  of  the  dif- 
ference equation.   Therefore  the  problem  of  numerical  sta- 
bility reduces  to  the  analysis  of  the  roots  of  the 
characteristic  equation  of  the  related  difference  equation 
of  (2-63). 
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C.   STABILITY  ANALYSIS  OF  A  ONE -STEP  METHOD 
Given  an  ODE 

y'  =  Ay  with   y(xQ)  =  yQ  (2-64) 

Assuming  that  A  is  constant,  the  analytic  solution  is  given 

by 

A(x  -x  ) 
Y(x)  =  yoe   n  °J 

Using  Euler  formula  (2-7),  the  related  difference  equation 
with  y^  =  Ayn  is 

y  Al  -  (1+Ah)y  =  0 
7n+l    v     J } n 

whose  characteristic  root  is 

u1  =  (1+Ah) 

Therefore  the  solution  is 

u*  =  (1+Ah)n 
Since  there  is  only  one  root,  (2-63)  reduces  to 

yn  =  d^  =  d^l+Ah)11 


where    dn    =   y    .      Therefore, 
1        J  o 


x    -x 
n      o 


Yn   =   yod+Ah)n   =   yQ(l+Ah)      n  (2-65) 

We  recall  from  calculus  that 

ex  =  lim   (l+h)x/h 
h+  0 
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and  using  the  fact  that  A  is  constant  (so  that  Ah  ->-  0  when 
h  ■*»  0)  ,  it  can  be  seen  that 

A(x  -x  )     A(x  -x  ) 
lin   (i+Ah)    n  °   =  e   n   °   . 
h  +  0  An 

It  follows  that  (2-65)  reduces  to 

y  =  y  e 

3  n   }  o 

The  method  is  then  (roughly  speaking)  stable,  which  means 
that  a  sequence  of  approximations  y  will  have  the  true 
solution  as  its  limit  if  the  corresponding  values  of  h  have 
zero  as  their  limit.   This  analysis  is  quite  simple  since 
no  extraneous  roots  represented  by  d7u?,  d,u,...,  are  intro- 
duced.  The  case  where  these  roots  are  introduced  will  be 
discussed  thoroughly  in  the  next  section.   The  only  error 
introduced  is  represented  by  the  particular  solution  of 
C2-59),  whose  behavior  has  been  thoroughly  studied  in 
Section  V-A  as  a  function  of  h. 

D.   STABILITY  ANALYSIS  OF  A  MULTI-STEP  METHOD 
Given  the  same  problem 

y*  =  Ay  with   y(xQ)  =  yQ 

the  corrector  equation  of  Milne's  P-C  set  (P-C-II)  is  used 
viz . 

y  o.i  =  y   i  +  \     y'  -.  +  4y'  +  y'     .  (2-66) 

7n+l    7n-l    3   7n-l    7n   7n+l  v     ' 

Since 

*A-1  =  Ayn-1 
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y '  =  Ay 


yn+i  =  A>Vi 

(2-66)  has  the  related  difference  equation  in  the  simplified 
form 


yn+l 


AT, 


Ah        4Ah 


J"  Xn  L- 


y 


n-l 


1  + 


Ah 


=  0. 


The  characteristic  equation  yields 


U- 


with  roots 


1  -  £l£i 


Ah 
3 


u- 


4Ah 


1  + 


Ah 


=  0. 


(2-67) 


ul  = 


2Ah  +/  9  +  3Ah' 
3  -  Ah 


Assuming  that  h  is  sufficiently  small,  the  binomial  theorem 
can  be  used  to  expand  the  square  root  and  then  simplifying, 
one  obtains 

ux  =  1  +  Ah  +  0(h) 
u2  =  -1  +  |  Ah  +  0(h) 


Since  n  =  (x  -x  )/h,  the  solution  of  (2-67)  can  be  written 


as 


(x  -x  ) 
(      -      ~)  /h  n  oJ 

yn=  d1(l+Ah+0(h))   n   °    +  d2(-l)n(l-  \   Ah+  0(h)) 


As  h  tends  to  zero,  this  solution  tends  to 


A(x  -x  )  -5-Afx  -x  ) 

,    vn   oy,j,nvn    3vn  oJ 
yn  =  d:e        +  d2(-l)   e 


(2-68) 
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Since  d,  is  given  by  the  initial  condition  and  d2  is  assumed 
to  have  been  supplied  by  an  auxiliary  one  step  method,  (2-68) 
can  be  completely  determined.   It  is  clear  that  there  is  an 
extraneous  root  introduced  as  a  result  of  the  use  of  a  2 
order  difference  equation  to  represent  a  first-order  dif- 
ferential equation.   The  root  u~  is  then  called  spurious, 
parasitic,  or  extraneous  and  has  no  relation  to  the   exact 
solution  of  the  differential  equation  but,  nevertheless,  is 
unavoidable.   This  is  clearly  seen  if  it  is  assumed  that 
d2  =  0,  in  which  case  (2-68)  reduces  to 

A(x  -x  ) 
v  n  oJ 
y   =  y   e 
7n   7  o 

This  is  indeed  the  true  solution  of  the  ODE.   But  in  general 
d~  f    0,  so  the  behavior  of  this  extraneous  root  remains  to 
be  studied,  since  it  affects  the  total  solution.   This  could 
be  done  proving  that,  if  A  >  0,  the  parasitic  solution 

u2  =  (-1)   e 

tneds  to  zero  exponentially  as  x   increases.   Then  the 
effect  of  the  spurious  solution  becomes  negligible  as  x 
increases  while  the  principal  solution 

A(x  -x  ) 
n      v  n  oJ 
U-.  =  e 

tends  toward  the  true  solution  values  yfx  ).   However  if 
A  <  0  the  spurious  solution  increases  exponentially  in  magni- 
tude and  alternates  in  sign  as  the  step-by-step  calculation 
progresses,  while  the  principal  root  u,  tends  to  one. 
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Moreover,  at  each  stage  of  the  calculation,  rounding  errors 
will  introduce  new  spurious  terms  of  the  same  type.   Since 
these  extraneous  roots  have  no  relation  to  the  exact  solu- 
tion, as  they  become  large,  they  will  dominate  the  results, 
thus  invalidating  the  numerical  solution.   This  situation  is 
referred  to  as  numerical  instability. 

The  same  analysis  can  be  generalized  by  going  back  to 
equation  (2-63) 

It  is  clear  then  that  the  principal  root  is  d-.u,  and 
d~u~,  d,u_ ,  .  .  .  d^+-.  u^.+  ,  are  extraneous  roots  that  are  intro- 
duced as  a  result  of  representing  a  first  order  differential 
equation  by  a  (K+l)    order  difference  equation.   By  analyz- 
ing the  principal  root  u,  which  has  the  form 

u1   =  1  +  Ah  +  0(h) 


or 


Ah     ,,  x 
u,  =  e    +  o(h) 


more  precise  requirements  about  numerical  stability  can  be 
determined.   Observe  that  u,  approximates  e     as  was  already 
shown. 

Comparing  (2-69)  to  the  exact  solution  of  y'  =  Ay,  viz, 

f      s     nAh 

yUn)  =  e 


it  is  clear  that  as  long  as  |u-,|  >  |u.|,  i  =  2,3, ...K+l, 


n 


when  n  increases,  the  extraneous  solution  represented  by  u . , 
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will  become  small  when  compared  to  the  principal  solution 
cLulj1.   In  this  situation  the  numerical  solution  is  expected 
to  be  stable.   However,  if  |u.|  >  | u, |  for  any  i  =  2,3,...K+1, 
u.  will  become  large  with  respect  to  u,  and  the  numerical 
solution  component  corresponding  to  u-  will  dominate  the 
solution.   This  situation  is  referred  to  as  numerical 
instability. 

Intuitively  then,  the  numerical  solution  will  be  valid 

if  I  u,  I  >  lu.l,  i  =  2,3,...K+1.   However,  recall  from  the 

1  1 '    '  1 '  '       *  '  ' 

previous  analysis  of  (2-59)  that  the  characteristic  equation 

(2-59)  for  the  accumulated  error  e  _,  ,  ,  is  exactly  the  same 

v     J  n+1 '  ' 

as  the  characteristic  equation  (2-63)  for  the  original 
difference  equation  of  (2-50).   For  a  valid  numerical  solu- 
tion, it  has  been  shown  that  the  total  error  e   must  not 

3  n 

grow  with  n.   From  (2-59),  it  was  observed  that 

enH  =  dlul  +  d2U2  +  •••  +  dK+lUK+l 

and  is  equivalent  to  the  condition  that  |u-|  <_   1,  i,2,3,...K+l 
It  is  then  possible  to  define  the  stability  of  a  method  in 
the  following  way. 

(i)   ABSOLUTELY  STABLE  if  |u.|  <  1,  i  =  1,2,...K+1 

(ii)   RELATIVELY  STABLE  if  lu.l  <  u, ,  i  =  2,3,...K+1 

K      J  '  l '  —   1 '       *  * 

Absolute  stability  does  not  imply  relative  stability.   In 
other  words,  a  numerical  solution  may  have  |u-|  £  1, 
i  =  1,2, ...K+l,  but  | u, |  <  |u.|  <  1,  i  ■  2,3, ...K+l. 
Summarizing  then,  it  can  be  seen  that  in  the  ODE 

y'  =  Ay 
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absolute  stability  and  relative  stability  are  both  important 
considerations  if  A  <  0  since  the  exact  solution  is  decreas- 
ing with  x  .   If  A  >  0,  the  exact  solution  is  growing  with 
x  ,  so  that  relative  stability  is  the  important  consideration 
In  other  words,  the  numerical  solution  is  valid  as  long  as  no 
component  of  u.  dominates  the  component  corresponding  to  the 
principal  root. 
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VI.   CONVERGENCE  OF  THE  CORRECTOR  IN  THE  P-C  METHODS 

The  term  convergence  has  been  mentioned  earlier  in  the 

study  of  the  propagation  of  the  error  as  a  function  of  h.   An 

analysis  of  the  convergence  properties  of  correctors  reveals 

that  this  convergence  indeed  depends  on  the  step  size  h.   A 

proper  selection  of  h  must  therefore  be  made  to  ensure  the 

convergence  of  the  corrector. 

th 


Let  y  +1  denote  the  j    approximation  of  the  soluti 


on 


value  at  x  _,,  .   Then  the  corrector  formula  can  be  written 
n+1 

in  the  form  with  $_,  f    0: 

y\t   =  y    +  h3  -,f(x  ^  ,y*L-, }  +  h   E  &.f(x   .  ,y   .) 
7n+l   7n-r     -1  k  n+l'^n+l7     *  =  o     n-i,/n-i-' 

(2-70) 

where  i  is  the  iteration  index.   Note  that  y  _,_-  =  y*\  ,  ,  and 
J  7n+l    7n+l' 

y-'J,  is  the  corrector  of  iterating  y \  ,  which  in  turn  is 
7n+l  6  7n+l 

the  corrector  of  (j-1)  iteration.   If  the  sequence  of  suc- 
cessive approximations  y  +1  ,  generated  by  iterating  the 
corrector,  converges  to  a  limit,  denoted  by  y  +, ,  then  yn+i 
satisfies  the  relation 

Cl  =  >Vr  +  B-lf(Vl>Cl'  +H  J0  6if(xn-i>>Vi) 

(2-71) 

Subtracting  (2-71)  from  (2-70)  yields 

yit\  -  Ci =  h6-i  [£txn+i>yn+i^  -  f'vrVi'1 

which  by  the  mean  value  theorem  can  be  written 
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yjtj  •  Cl  =  hB-l  [yn+l  "  d]  fy(Xn+l'7n+l} 

*        i+1      ii 
where  yn+1  lies  between  yn+1  and  yn+n«   If  I f v I  £  C,  where 

C  is  a  positive  constant,  in  the  neighborhood  of  (x  .  n  ,y  _,,  ) 
r  7  &  v  n+1   n+ly 

then 

i+1     *   ,     ■    iii       * 
vJ-v      <hfi    C  vJ    -v 
|yn+l   yn+l!  -  n|P-l|L|yn+l   yn+l'' 

It  follows  by  induction  that 

|yj+1  -  y*   I  <  rhlB   IClj+1  |y°    -  y*   I 
|yn+l   yn+ll  -  inlp-llLJ     lyn+l   yn+l! 

Therefore,  the  iteration  of  the  corrector  will  converge  to 

* 

the    limit   y    , ,    provided   that 
7n+l   r 

|h.e_1C|    <    1  (2-72) 

With  this  expression  the  conditions  for  convergence  of  the 

various  corrector  formulas  considered  can  be  written  down 

immediately: 

Condition  for 


Corrector 

1/2 

Com 
hC 

.L1UII    J-\J   X 

mergence 

Euler 

<  2 

Milne 

1/3 

hC 

<  3 

Hamming 

3/8 

hC 

<  8/3 

Third-order  Adams 

5/12 

hC 

<  12/5 

Fourth-order  Adams 

3/8 

hC 

<  8/3 

Both  the  Nystrom  and  second  order  Adams  P-C  sets  used  the 
Euler  corrector  while  Hermite  used  the  Milne  corrector  and 
thus  are  included  in  the  above  expressions. 
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An  important  point  that  must  be  brought  out  here  is  that 
convergence  will  only  assure  that  the  sequence  of  iterations 

y ■ . n  will  converge  to  some  definite  value  y  .  n ,  but  not 

' n+1  to  J  n+1 

necessarily  to  the  true  solution  y (x  ).   Stability  is  still 
the  yardstick  of  a  valid  numerical  solution   though  (2-72) 
provides  a  good  rule  of  thumb  to  follow  in  the  proper  selec- 
tion of  h.   Indeed,  later  it  will  be  shown  experimentally 
that  absolute  stability  can  not  be  attained  beyond  the  bound 
for  h  given  by  (2-72)  since  convergence  is  a  necessary  but 
not  a  sufficient  condition  for  stability  of  a  numerical 
solution  (cf.  Ralston  [Ref.  6]). 
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VII.   INFLUENCE  OF  THE  PROPAGATED  ERROR 

In  order  to  test  the  idea  presented  in  the  previous  sec- 
tion dealing  with  the  propagation  of  error,  the  numerical 

example  given  in  Section  V-A  was  run  on  an  IBM  360/67  compu- 

-  7 
ter  where  the  round-off  err  bound  is  10   using  single 

-14 
precision  (SP)  and  10     in  double  precision  (DP) .   Recall 

that  the  ODE  used  was 

y1  =  -y  with  y(0)  =  1 

and  the  final  error  expression  using  the  Euler  predictor 
formula  at  x  =  1  is 

A   series  of  h  values  from  h  =  10    to  h  =  1.0  was  used.   The 
numerical  values  of  e   were  computed  both  in  single  precision 

and  double  precision.   Table  1  shows  the  actual  numerical 

_  3 
data  obtained.   These  data  show  that  for  large  h  to  h  ^  10   , 

the  SP  and  DP  results  are  essentially  identical,  but  as  h 

decreases  the  round-off  errors  tend  to  increase  while  the 

truncation  error  tends  to  zero.   The  effect  of  the  round-off 

error  in  DP  remains  negligible  since,  in  double  precision, 

14  digit  accuracy  is  obtained.   These  results  confirm  the 

idea  that  the  errors  decrease  with  decreasing  values  of  h 

to  a  certain  point,  beyond  which  the  errors  increase.   A 

plot  of  the  data  obtained  is  shown  in  Figure  2. 
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TABLE  1 

INFLUENCE  OF  PROPAGATED  ERROR  (Z) 
FOR  y'  =  -y  AT  x=1.0 


h 

7* 
^ESP 

EEDP 

l'lO"6 

1.595-10"1 

8.763-10"7 

1-10"5 

1.596-10"2 

8.593-10"6 

5-10"5 

3.434-10"3 

4.295-10"5 

-4 
1-10 

1.780-10"3 

8.590-10"5 

5-10"4 

9.011-10"4 

2.147-10"4 

-4 
5*10    4 

7.722-10"4 

4.294-10"4 

1-10"3 

1.029-10"3 
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ESP  -  Error  in  Single  Precision 
EDP  -  Error  in  Double  Precision 


Figure  2.   Influence  of  Propagated  Error. 
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VIII.   STABILITY  BOUNDS  OF  P-C  SETS 

Prior  discussion  about  the  stability  of  the  predictor- 
corrector  methods  relates  only  to  the  limiting  properties 
of  the  roots  of  the  characteristic  equation  as  the  interval 
of  integration  approaches  zero.   In  many  applications, 
however,  additional  information  about  the  actual  size  of  h 
is  needed  to  ensure  stability  of  a  method.   The  number  of 
results  discussed  in  published  papers  on  this  subject  is 
extensive.   But  the  interesting  point  observed  is  that 
although  they  all  started  with  the  analysis  of  the  roots  of 
the  characteristic  equations,  they  varied  quite  extensively 
in  the  mode  of  the  predictor-corrector  applications  and  in 
the  number  of  P-C  sets  considered.   To  be  more  specific, 
some  published  papers  are  cited.   The  real  negative  stability 
expression  used  was  hX  where  X  =    f   in  a  single  ODE.   Chase 
[Ref.  7]  analyzed  the  Milne  P-C  mode  without  iteration.   The 
The  real  negative  stability  bounds  found  were  -0  .  8 <  hX  <  -  0 . 3  , 
while  on  the  other  hand  using  the  same  P-C  mode,  but  now 
iterated  to  convergence,  resulted  in  numerical  instabilities 
for  all  negative  X.   Hamming  [Ref.  4]  used  three  modes  of 
his  P-C  set.   First,  iterated  to  convergence,  the  resulting 
stability  bounds  were  -0.5  <  hX  <  0;  second,  with  truncation 
error  modification  but  without  iteration  the  results  were 
-0.85  <  hX  <  0;  third,  still  with  truncation  error  modifica- 
tion but  with  only  two  iterations,  the  results  were 
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-0.9  <  hX    <    0.   Lapidus  and  Seinfeld  [Ref.  2]  using  the  PECE 
and  PE(CE)   methods,  where  s  is  the  iteration  index,  pub- 
lished the  real  negative  stability  bounds  for  all  available 
predictor-corrector  methods.   Their  algorithm  is  quite 
different  from  the  algorithm  used  in  this  paper  in  that  a 
final  derivative  evaluation  is  computed  before  terminating 
the  computation  with  or  without  iteration.   The  symbol  PECE 
or  PE(CE)   is  used  to  denote  their  algorithms.   The  algorithm 
considered  in  this  paper  does  not  require  a  final  derivative 
evaluation  before  termination  with  or  without  iteration,  and 
hence  can  be  denoted  by  P(EC)  . 

To  establish  the  real  negative  stability  bounds  for 
the  P-C  sets  considered  herein,  the  experimental  procedure 
used  by  Chase  to  obtain  the  real  negative  stability  bounds 
for  hX   will  be  followed.   Chase  used  the  test  ODE 

y'  =  100  -  lOOy  with  y(0)  =  0 

Using  different  values  of  h,  he  analyzed  the  behavior  of  the 
error  until  actual  instability  occurred.   The  P-C  mode  will 

be  different  in  the  sense  that  a  standard  application  involv- 

2 
ing  two  iterations,  P(EC)  ,  will  be  used  for  all  numerical 

methods  considered.   The  choice  of  just  two  iterations  was 

incluenced  first  by  the  published  paper  of  Hull  and  Creemer 

[Ref.  8].   They  conducted  an  experiment  using  the  P(EC) 

mode  applied  to  the  Adams'  methods  only.   After  analyzing 

several  numerical  results,  they  concluded  that  the  best 

method  is  S  =  2 ,  based  on  the  cost  of  computation,  average 

accuracy,  and  stability.   Second,  the  choice  of  iterations 
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was  influenced  by  the  analysis  of  the  corrector  in  terms  of 
minimum  truncation  error.   Consider  Euler's  P-C  set  (P-C-I), 
which  has  the  form 

Predictor:   y  .-  =  y   +  hy ' 
7n+l    ;n  7n 

Corrector:   yR+1  =  yn  +  \   [y^+1  +  y^] , 

as  applied  to  the  ODE  y'  =  Ay. 

Ignoring  the  corrector,  the  single  characteristic  root 
is  obtained  from  the  predictor  as  follows: 

>Vl  -  >^n  -  h^n  =  ° 
u1    -  (1  +  hA)  =  0 

u,  =  1  +  hX. 

Now  apply  the  corrector  once.   Substituting  the  predictor 
into  the  corrector  yields 

2 


=  (l  +  hA  +  Uf-)    yn 


y 


The  single  characteristic  root  is 

u,  ■  1  +  hA  +  iMii 

1  2 

Now  apply  the  corrector  twice  obtaining 

2  2 

yn+i  ■  rn  + 1 1^1  +  hx  +  V-)  yn  +  XV 

The  characteristic  root  is 
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n   .,   h2X2   h3X3 
ul  =  1   +  hX  +  ~2~   +  T- 

Continuing  this  process  yields  these  results 

0  corrector:   u,  =  1  +  hX 

h2X2 

1  corrector:   u,  =  1  +  hX  +  —- — 

2  correctors :  u,  =  1  +  hX  +  — s —  +  — ■? — 

2    2  3    3  4    4 

■7  4.  iu-vhA  h   X  hX 

3  correctors:    un    =    1   +   hX   +  — = —  +   — : —  +  —5 — 

1  Z  4  o 

However,    the    true    solution   of  y1    =    X      is    given   by 

eh*    =    1    +    hX    +   4^~  +   4^  +   44  +    0(h5) 

Z  o  Z4 

and  thus  the  truncation  error   T(x,h) ,  given  by 

hX 
u,  -  e 

is,    first  with    0    corrector 

2    2 
T(x,h)    =    (1  +  hX)    -    (1+hX  +  ^~—  +    0(h3)) 

T(x,h)    =   ^*         -    0(h3) 


By  continuing  this  procedure  with  the  1,2,3  correctors  and 
ignoring  higher  order  terms,  there  results 

0  corrector:  T(x,h)  =  -X2h2/2 

1  corrector:  T(x,h)  =  -X3h3/6 

2  correctors:  T(x,h)  =  X3h3/12 

3  correctors:  T(x,h)  =  X3h3/12 
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It  can  be  seen  that  the  use  of  the  corrector  more  than  twice 
seems  to  yield  little  advantage  in  this  case.   Thus  the  mode 

of  P-C  application  to  be  used  will  be  represented  by  the 

2 
symbol  P (EC)  .   In  this  mode  of  application  the  stability  of 

the  P-C  set  depends  on  both  the  predictor  and  corrector 
applications,  though  more  so  on  the  corrector  equation.   This 
assertion  will  be  clearly  illustrated  in  the  case  of  the 
Milne,  Hamming  and  Nystrom  P-C  sets.   Thus,  a  priori  know- 
ledge of  the  stability  of  a  method  can  be  obtained  by 
analyzing  the  roots  of  the  characteristic  equation  of  the 
corrector  formula. 

A.   EXPERIMENTAL  PROCEDURE 

Two  general  algorithms  were  developed  to  implement  the 
predictor-corrector  methods  in  Section  IV-E.   The  first 
algorithm  covers  seven  of  the  P-C  sets,  while  the  second 
algorithm  covers  the  8    P-C  set,  i.e.,  the  Hamming  modified 
predictor-corrector  set.   The  flow  chart  for  ALGORITHM  1  is 
on  page  164  and  that  for  ALGORITHM  2  is  on  page  167.   How- 
ever, the  general  steps  for  each  algorithm  will  be  outlined 
here . 

Given  the  necessary  starting  values,  the  step  size,  the 
convergence  term  and  the  range  of  integration,  the  iterative 
predictor-corrector  method  for  computing  the  numerical  solu- 
tion y_+]  is  set  up  as  follows: 

ALGORITHM  1: 

Step  1:   Compute  the  predicted  value  by  the  predictor 
formula. 
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Step  2:   Compute  the  derivative  of  the  predicted  value. 

Step  3:   Compute  the  corrected  value  by  the  corrector 
formula. 

Step  4:   Test  for  convergence.   If  convergence  is  attained 
go  to  step  7,  otherwise  proceed  to  the  next  step. 

Step  5:   Compute  the  derivative  of  the  corrected  value. 

Step  6:   Compute  the  new  corrected  value  by  the  corrector 
formula,  using  the  value  computed  in  step  3  as 
the  new  predicted  value. 

Step  7:   The  result  is  taken  as  the  desired  solution 

point.   Advance  the  integration  point  by  the  step 
size.   Return  to  step  1  to  compute  the  next 
solution  point. 


ALGORITHM  2 


Step  1:   Compute  the  predicted  value  by  the  predictor 
formula. 

Step  2:   Modify  the  predicted  value  by  adding  the  trun- 
cation error*  of  the  predictor  formula. 

Step  3:   Compute  the  derivative  of  the  modified  predicted 
value . 

Step  4:   Compute  the  corrected  value  by  the  corrector 
formula. 

Step  5:   Modify  the  corrected  value  by  adding  the  trun- 
cation error*  of  the  corrector  formula. 

Step  6:   Test  for  convergence.   If  convergence  is  attained 
go  to  step  10;  otherwise  proceed  to  the  next 
step. 

Step  7:   Compute  the  derivative  of  the  modified  corrected 
value . 

Step  8:   Compute  the  new  corrected  value  by  the  corrector 
formula  using  the  value  as  computed  in  step  5 
as  the  new  modified  predicted  value. 


* 

Truncation  error  as  used  in  the  formula  is  not  the 
original  truncation  error  but  the  modified  truncation  term 
expressed  as  a  function  of  the  difference  between  the  cor- 
rected and  predicted  values. 
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Step  9:   Modify  the  corrected  value. 

Step  10:  The  result  obtained  is  accepted  as  the  desired 
solution  point.  Advance  the  integration  point 
by  the  step  size.  Return  to  step  1  to  compute 
the  next  solution  point. 

2 
In  both  algorithms  the  mode  of  application  is  P(EC)   as  can 

be  readily  verified.   In  this  scheme,  if  the  corrector  is 

used  only  once,  then  the  number  of  function  evaluations  needed 

is  also  one,  since  it  is  assumed  that  the  function  value  for 

the  corrector  has  already  been  computed;  if  the  corrector  is 

used  twice  then  two  function  evaluations  are  needed.   Thus 

the  number  of  function  evaluations  is  determined  by  the 

number  of  iterations. 

The  P-C  algorithms  were  coded  using  the  FORTRAN  language. 
The  fortran  programs  for  P-C-I  to  P-C-VIII  are  listed  on 
pages  170  to  200.   In  each  program  the  meanings  of  symbols 
and  parameters  are  explained  through  narrative  comments. 

To  determine  whether  the  stability  limits  have  been 
reached,  certain  criteria  must  be  followed.   Often,  when  the 
stability  bound  is  approached  through  increasing  step  size  h, 
the  method  begins  to  lose  accuracy.   An  inaccurate,  though 
stable,  solution  may  often  result  in  oscillations  which 
appear  at  first  to  be  due  to  actual  numerical  instability. 
Actual  instability,  associated  with  too  large  a  value  of  h, 
usually  results  in  an  increasing  oscillation  in  the  error, 
or  simply,  the  growth  of  the  error  in  one  direction.   The 
particular  error  behavior  depends  on  the  sign  of  the  para- 
sitic root  causing  the  instability;  a  negative  root  causes 
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oscillations  and  a  positive  root  causes  a  uniform  error 
growth. 

B.   PREDICTED  STABILITY  CHARACTERISTIC  OF  THE  P-C  SETS 

To  gain  a  priori  knowledge  of  the  behavior  of  the  sta- 
bility characteristics  of  the  P-C  sets,  the  roots  of  the 
characteristic  equations  of  the  corrector  formulas  will  be 
studied.   It  must  be  noted  that  the  behavior  of  the  charac- 
teristic roots,  though  dominating  the  P-C  method  used  as  a 

set,  will  be  influenced  by  the  behavior  of  the  predictor 

2 
formula.   Since  the  chosen  mode  of  application,  P(EC)  ,  is 

quite  different  from  all  the  other  modes  of  applications 
published  in  different  papers,  the  stability  limits  obtained 
by  several  authors  could  not  be  used  as  the  stability  limits 
of  the  P-C  set  considered  here.   Thus  the  actual  real  nega- 
tive stability  limits  will  be  determined  through  numerical 
experiment.   However,  as  noted  previously,  the  analysis  of 
the  characteristic  roots  remains  the  same,  and  as  such  the 
published  results  are  used  and  references  are  cited.   This 
is  done  solely  to  conserve  space  and  to  avoid  the  tedious, 
repetitious  process  needed  to  solve  the  characteristic 
equations.   The  procedures  have  been  amply  presented  and 
illustrated  in  detail  in  the  section  on  numerical  stability 
of  predictor-corrector  methods. 

a)   Analysis  of  the  characteristic  roots  of  the  corrector 
formulas . 

The  P-C  sets  with  common  correctors  will  be  grouped 
together  except  for  the  case  of  the  second  order  Adams' 
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method  which  has  the  same  corrector  as  Euler,  and  which  will 
be  included  in  Adams'  type. 

Euler  (P-C-I)  and  Nystrom  (P-C-III): 
The  corrector  formula  is  given  by 


Vl  =  yn  +  2    [yn+l  +  K] 
The  related  difference  equation  is  in  the  form 

(1  '  ¥^  yn+l  "  (1  +  ¥}  yn  =  °«  Where  C  =  VX,7) 

The  single  characteristic  root  is 

hC 


1  + 


2 


Ul  = 


1  - 


hC 


It  can  be  seen  that  if  C  <  0,  (i.e.,  the  derivative  3f/8y  is 
negative),  then  for  absolute  stability  to  occur,  |u, |  <  1 
which  is  satisfied  if  hC/2  <  1.   Thus  the  numerical  solution 
will  decrease  as  does  the  exact  solution   If  C  >  0,  as  long 
as  hC/2  <  1  then  the  numerical  solution  increases  as  does  the 
exact  solution.   Figure  3  exhibits  the  behavior  of  the  root 
versus  hC  as  shown  by  Lapidus  and  Seinfeld  [Ref.  2]. 

Milne  (P-C-II)  and  Hermite  (P-C-IV): 
The  corrector  is  of  the  form 


yn*i =  >vi +  3  [yA-i +  4^A +  W 


The  related  difference  equation  is 


[• 


hC 


yn+l 


4hC 


y 


n 


1  + 


hC 


7n-l 
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hC 

Figure  3.   Characteristic  Root  for  Euler  Corrector  Formula 
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The  roots  of  the  characteristic  equation,  as  derived  in 
section  V-D,  are  of  the  form 

u1   =  1  +  hC  +  0(h) 

u2  =  -1  +  i  hC  +  0(h) 

which  can  be  written  as 

hC 
u-.  -  e 


u2  -  e 


hC/3 


If  C  >  0,  u1  behaves  like  the  exact  solution  and  u~  dies  out 
since,  |u2|  <  1.   When,  however,  C  <  0 ,  u,  decreases  as 
does  the  exact  solution,  but  u~  increases.   Thus  the  corrector 
is  relatively  stable  but  for  C  >  0,  has  no  real  negative 
stability  bounds.   Figure  4  shows  the  behavior  of  the  roots 
as  given  by  Chase  [Ref .  7] . 

Hamming  (P-C-V) : 

The  corrector  is  given  by 

The  related  difference  equation  is  of  the  form 


1)7n+1    +    l 


3hC         ,1„  .    |9    +    3hC  .       3hC  .1  =    0 

8  4      yn  8      yn-l        8   yn-2 


The  characteristic  equation  is 


u3   3hC  .    +  u2   9  +  3hC   .    3hC    1  .  „ 


Chase  [Ref.  7]  analyzed  the  root  behavior  of  this  character 
istic  equation.   Figure  5  shows  the  behavior  of  these 

75 


1.6 


1.4 


1.2  <- 


1.0 


0.8  - 


0.6  — 


0.4  '- 


0.2  '- 


Figure  4.   Characteristic  Roots  for  Milne  Corrector  Formula 
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1.  D 

x  -  Positive  Roots 
N  -  Negative  Roots 
0  -  Complex  Roots 
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— * 

n                                        ^ 
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1        1 

1           1           1 

-2.4    -1.6    -0.8  hc   0     0.8     1.6     2.4 
Figure  5.   Characteristic  Roots  for  Hamming  Corrector  Formula 
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characteristic  roots  versus  hC.   Analyzing  the  graph  of 
Figure  5,  it  is  evident  that  the  corrector  is  stable  for 
-2.4  < hC  <  0.   This  corrector  has  also  a  wide  range  of 
relative  stability  as  Chase  has  shown.   In  the  graph  the 
symbol  used  for  positive  roots  is  an  x,  that  for  negative 
roots  is  an  N,  and  for  the  complex  roots  is  an  0. 

ADAMS-MOULTON  CORRECTORS  (Second  Order  (P-C-VI,  Third 
Order  (P-C-VII) ,  Fourth  Order  (P-C-VIII)): 

Since  these  formulas  are  all  of  the  same  type,  they  will 

be  studied  together  and  extended  to  the  general  case  of  the 

Adams-Moulton  types.   For  the  second  order,  the  corrector  is 

given  by 

>Vi =  ^  + 1  [K+i+  yk]- 

The    related   difference    equation    is    of   the    form 

[1    "   1^]    y   a.i    '    [1   +   t^I    y      =0 
1  2    J    7n+l         l  2        7n 

with  characteristic  equation 

[1  -  |£]  u  -  [1  +  |£]  =  0 

In  the  limit,  as  h  -*■  0 ,  the  root  becomes 


u,  =  1 . 


For  the  third  order,  the  corrector  is 


yn+i  -  yn  + 12  i^i +  8^A  -  y;-i] 


The  related  difference  equation  is  given  by 
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t1  -  fi hcI  >vi  -  t1  +  !hc)  yn  +  t§  vi  -  ° 

with  characteristic  equation 

[1  -  jj   hC]  u2  -  [1  +  |  hC]  u  +  yj  =  0 

In  the  limiting  case  h  -*•  0 ,  the  roots  are 

u(u-l)  =  0 

u,  =  1  ;  u~  =  0  . 

For  the  fourth  order,  the  corrector  is 

The  related  difference  equation  is 

t1  "  Hr>  ^n+l  "  H  +  IT  hC]^n  +  if  Xn-1  #  y  _,-  0 


24   7n-2 


with  characteristic  equation 


[1  -  |  hC]u3  -  [1  +  i|  hC]u2  +  |^  u  -  ^  =  0 

In  the  limiting  case  as  h  ->  0,  the  roots  are 
u2(u-l)  =  0 
u,  =  1,   u~  =  u.,  =  0 

Thus,  it  can  be  seen  that,  in  each  case,  the  dominant  root 
is  1  while  all  the  other  roots  (the  parasitic  ones)  lie  at 
the  origin  when  h  =  0.   In  the  general  q-step  Adams -Moulton 
formula  the  equivalent  characteristic  equation  is 

uq_1(u-l)  =  0  . 
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Thus,  in  every  case,  the  A-M  formulas  have  one  root  on  the 
unit  circle  with  all  others  at  the  origin  when  h  =  0. 
As  such  they  would  seem  to  have  excellent  stability  charac- 
teristics no  matter  what  q-step  formula  is  used.   When  h 
increases  from  zero,  the  parasitic  roots  move  toward  the 
unit  circle.   However,  significantly  large  values  of  h  can 
be  reached  before  instability  occurs.   Crane  and  Klopfen- 
stein  [Ref.  9]  analyzed  the  behavior  of  the  general  Adams- 
Moulton  formula  and  showed  that  it  has  stability  bounds 
-1.3  <  hC  <  0 .   It  is  interesting  to  note  that,  in  general, 
as  the  order  of  the  Adams  corrector  increases,  accuracy 
increases  but  stability  decreases.   This  is  the  main  reason 
why  higher  order  formulas  are  not  considered  (i.e.,  fifth, 
sixth,  seventh,  eighth).   The  graph  of  the  characteristic 
roots  of  the  general  A-M  formulas,  as  presented  by  Crane  and 
Klopfstein,  was  not  reproduced  since  it  would  give  no 
additional  information. 

In  the  foregoing  analysis  of  the  characteristic  roots, 
the  predicted  stability  limits  served  as  a  guide  in  the 
choice  of  the  step  size  h  to  be  used  in  the  numerical 
experiment  for  obtaining  the  actual  real  negative  stability 
bounds  for  each  P-C  set  considered 

C.   NUMERICAL  RESULTS  FOR  REAL  NEGATIVE  STABILITY  LIMITS 

In  order  to  obtain  actual  real  negative  stability  bounds 
of  the  P-C  sets  considered  in  the  P(EC)   mode  of  application, 


the  ODE 


was  chosen,  where  fy(x,y)  =  C  =  -1  <  0. 


y'  =  -y  with  y(0)  =  1 
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Each  P-C  set  was  run  with  increasing  values  of  step 
size  h,  in  small  increments,  to  determine  more  precisely  the 
actual  stability  limits.   The  criteria  for  actual  numerical 
instability  set  forth  previously,  were  followed  in  the 
analysis  of  data  obtained.   The  prior  knowledge  of  the  roots 
of  the  characteristic  equation  obtained  in  the  previous 
analysis  narrowed  the  choice  of  starting  values  for  the  step 
size  h  for  each  method.   Each  method  was  run  in  single 
precision  because  of  the  large  values  of  h. 

1.   Euler  (P-C-I) 

Since  the  predicted  stability  of  the  corrector  is 
in  the  range  -2  <  hC  <  0,  the  starting  value  chosen  for  h 
was  1.1,  then  was  increased  by  0.1  until  instability  occurred 
Table  2  shows  selected  results  for  different  values  of  h. 
At  h  =  1.1,  it  is  absolutely  stable  since  the  error  continues 
to  decrease  as  x  increases.   The  first  sign  of  oscilla- 
tion was  observed  at  h  =  1.3  but  it  is  still  stable  because 
the  error  is  decreasing.   Oscillation  becomes  evident  at 
h  =  1.6,  but  although  the  solution  is  inaccurate,  it  is 
not  yet  unstable,  because  there  is  no  uniformly  increasing 
trend  of  oscillation.   Actual  numerical  instability  occurred 
at  h  =  1.9  as  is  clearly  shown  by  the  uniformly  increasing 
oscillation.   Thus  this  method  is  stable  within  -1.9<hC<  0. 
As  compared  to  the  predicted  stability  region  -2.0  <  hC  <  0 
obtained  in  [Ref.  2]  where  in  the  corrector  is  iterated  to 
covergence,  this  was  expected  because  of  the  use  of  only 
two  iterations. 


TABLE  2 

ABSOLUTE  ERROR  (Z )  FOR  EULER  (P-C-I) 

ODE:   y'  =  -y  with  y(0)  =  1 
Single  Precision 


h  =  1.1 

h  =  1.3 

h 

x              E 

1.1        -2.239-10"2 

1.3         -8.023-10"2 

2.2        -5.765-10"2 

2.6         -1.019-10"1 

3.3         -2.321'10~2 

3.9         -2.748*10"2 

4.4        -1.610'10~2 

5.2         -2.957*10"2 

5.5        -6.078'10"3 

6.5        -1.816*10"3 

6.6        -3.421-10"3 

7.8        -8.301*10"3 

7.7        -1.266-10"3 

9.1         1.629'10"3 

8.8        -6.549-10"4 

(STABLE) 

9.9        -2.406*10"4 

(STABLE) 

h  =  1.6 

h  =  1.9 

x              E 

x              E 

1.6         -2.733-10"1 

1.9        -6.696-10"1 

3.2         -6.44-10"2 

3.8          5.308-10"1 

4.8         -1.449-10"1 

5.7         -1.546 

6.4         7.894-10"2 

7.6          3.020 

8.0        -1.442-10"1 

9.5         -6.363 

9.6         1.634-10"1 

(UNSTABLE) 

(STABLE  BUT  INACCURATE) 
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2.  Milne  (P-C-II) 

A  priori  knowledge  indicates  that  this  method  has  no 
real  negative  stability,  so  the  starting  value  of  h  was 
chosen  as  0.1.   Table  3  shows  the  selected  data  obtained. 
Even  at  the  starting  value  at  h  =  0.1,  the  buildup  of  error 
is  consistent  though  small.   Then  at  h  =  0.4  oscillations 
become  evident  and  the  error  tends  to  increase  as  x  becomes 
large.   At  h  =  0.7  the  same  increasing  error  tendency  was 
observed.   These  observations,  based  on  actual  results, 
confirmed  the  predicted  instability  of  this  method  for  C  <0. 
Thus  it  can  be  concluded  that  this  method  has  no  real 
negative  stability  bounds,  as  is  evident  from  the  actual 
data  presented  on  Table  3. 

3.  Nystrom  (P-C-III) 

The  selected  starting  value  is  close  to  the  value 
used  in  Euler  since  this  method  has  the  same  corrector.   The 
program  for  Nystrom  was  run  starting  at  h  =  1.0  in  increments 
of  0.1  until  instability  occurred.   Table  4  shows  the  actual 
data  obtained.   The  result  for  h  =  1.1  to  1.3  show  absolute 
stability.   For  h  =  1.5,  the  error  starts  oscillating  but  not 
at  a  uniformly  increasing  rate.   The  solution  corresponding 
to  this  data  is  inaccurate  but  stable.   This  inaccuracy 
continues  up  to  h  =  1.7.   For  h  =  1.8  actual  numerical 
instability  has  occurred  as  shown  by  the  increasing  oscilla- 
tion.  Thus  this  method  is  stable  for  -1.8  <  hC  <  0 .   This 
range  of  stability  is  less  than  that  for  the  Euler  method 
though  they  have  the  same  corrector  equation.   This  confirmed 
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TABLE  3 

ABSOLUTE  ERROR  (E)  FOR  MILNE  (P-C-II) 

ODE:   y'=-y  with  y(0)=l 
Single  Precision 


h  =  0.1 

h  =  0.4 

x              M 

x              M 

1.0           1*10"7 
2.0           5-10"7 
3.0           5-10"7 
4.0           5*10"7 
5.0           6-10"7 
6.0           1-10"6 
7.0         1.6-10"6 
8.0         2.1-10"6 
9.0         3.0-10"6 
10.0         4.2*10"6 
(UNSTABLE) 

2.4         5  . 1 7  '  1 0  ~  5 
3.6        -5.86-10"5 
4.8          1.115'10~4 
6.0         -1.555-10"4 
7.2          2.393-10"4 
8.4         -3.58-10"4 
9.6          5.394*10"4 
(UNSTABLE) 

h  =  0.7 

x             EM 

2.8        2.584-10"4 

4.2        7.106-10"4 

5.6         8.312-10"4 

7.0        1.186-10"3 

8.4         1.918-10"3 

9.8         3.213*10"3 
(UNSTABLE) 
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TABLE  4 
ABSOLUTE  ERROR  (£)  FOR  NYSTROM  (P-C-III) 


ODE 


y 


y  with  y(0)=l 


Single  Precision 


h  =  1.1 

h  =  1.3 

x             EN 

E  XT 

x              N 

2.2  3.465-10"2 

3.3  4.369-10"2 

4.4  3.043-10"2 

5.5  1.724-10"2 

6.6  7.181*10"3 

7.7  3.062-10"3 

8.8  4.488-10"4 

9.9  3.363-10"4 
(STABLE) 

2.6        5.570-10"2 
3.9         5.976*10"2 
5.2         5.013'10"2 
6.5        1.446*10"2 
7.8        1.914*10"2 
9.1       -1.183*10"2 
(STABLE) 

h  =  1.5 

h  =  1.8 

x              N 

^m 

x              N 

3.0         7.005-10"2 
4.5         6.538'KT2 
6.0         8.581*10"2 
7.5        -2.784'10"2 
9.0         1.128-10"1 
(STABLE  BUT  INACCURATE) 

3.6       -1.865-10"2 
5.4         2.057"10"1 
7.2        -3.134'10"1 
9.0        1.132 
10.8       -3.2066 
(UNSTABLE) 
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the  idea  that  the  stability  of  the  predictor-corrector 
depends  on  both  the  predictor  and  corrector  formula. 

4.  Hermite  (P-C-IV) 

This  method  has  the  same  corrector  as  the  Milne 
method,  hence  the  choice  of  starting  value  is  the  same, 
h  =  0.1,  which  is  then  Incremented  by  0.1  until  instability 
occurred.   Table  5  shows  the  selected  actual  data  obtained. 
From  these  results  it  is  evident  that  this  method  is  unstable 
since  for  all  values  of  h  the  error  growth  is  steadily 
increasing.   This  is  expected  since  it  was  predicted  that 
the  corrector  has  a  dominant  role  in  the  stability  of  this 
method  as  a  P-C  set.   Thus  this  method,  which  uses  a  Milne 
corrector,  followed  the  instability  behavior  of  the  Milne 
corrector,  within  two  iterations.   Therefore,  it  can  be 
seen  that  the  Hermite  P-C  set  has  no  real  negative  stability 
bounds . 

5.  Hamming  CP-C-V) 

The  chosen  starting  value  was  h  =  0.5  with  increment 
0.1.   The  predicted  stability  of  the  corrector  is  -2.4  < hC <  0 
Table  6  shows,  that  for  h  =  0.5  to  h  =  0.6,  the  method  is 
absolutely  stable.   For  h  =  0.7  to  h  =  0.8  the  method  is 
stable  although  the  numerical  solution  becomes  inaccurate. 
For  h  =  0.9  the  error  grows  in  one  direction  thus  resulting 
in  actual  instability.   Thus  the  method  is  stable  within 
the  range  of  -0.9  <  hC  <  0.   Compared  to  the  predicted  sta- 
bility of  the  Hamming  corrector,  given  by  -2.4  <  hC  <  0,  the 
range  of  stability  is  greatly  reduced.   This  is  due  to  the 
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TABLE  5 

ABSOLUTE  ERROR  (E)  FOR  HERMITE  (P-C-IV) 

ODE:   y'=-y  with  y(0)=l 
Single  Precision 


:: — t:  -,  v : '  ' ^— : : ,[  ~  : — — : — — 

h  =  0.2 

■— — T-" '    ~  ■ ■ ,'-   ■  ■  T  , ■  ,  , 

h  =  0.5 

x             EHE 

x             EHE 

1.0         -l'lO"6 
2.0          6*10"6 
3.0          -5.1-10"6 
4.0          9.1-10"6 
5.0         -1.18*10"5 
6.0          1.71-10"5 
7.0         -2.4-10"5 
8.0          3.39*10~5 
9.0         -4.78*10"5 
10.0          6.74-10"5 
(UNSTABLE) 

1.0        2.820-10"4 
2.0        3.966-10"4 
3.0        6.146-10"4 
4.0        9.347-10"4 
5.0        1.401'10"3 
6.0        2.088-10"3 
7.0        3.106'10"3 
8.0        4.617-10"3 
9.0        6.863'10"3 
10.0        1.019*10"2 
(UNSTABLE) 

h  =  0.8 

x             EHE 

1.6         2.388'10"3 
3.2         6.331-10"3 
4.8        1.521*10"2 

6.4         3.298*10"2 

8.0         7.033U0"2 

9.6        1.496-10"1 
(UNSTABLE) 
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TABLE  6 

ABSOLUTE  ERROR  (I)  FOR  HAMMING  (P-C-V) 

ODE:   y'=-y  with  y(0)=l 
Single  Precision 


h  =  0.5 

h  =  0.6 

x             ZHA 

x             ZHA 

2.0        -2.056-10"4 
3.0        -1.418-10"4 
4.0        -8.06'10"5 
5.0        -4.15*10"5 
6.0        -2.02-10"5 
7.0        -0.6-10"6 
8.0        -4.6*10"6 
9.0        -2.2*10"6 

10.0       -i.no-6 

(STABLE) 

2.4         -3.741*10"4 
3.6         -2 . 897  * 10"4 
4.8        -1.689-10"4 
6.0        -9.37-10"5 
7.2         -5.62'10"5 
8.4        -3.78*10"5 
9.6        -2.81-10"5 
(STABLE) 

h  =  0.7 

h  =  0.8 

h  =  0.9 

x        EHA 

x        ZHA 

x        ZHA 

2.8    -5.709-10"4 
4.2    -5.717*10~4 
5.6    -3.944-10"4 
7.0    -3.047"10"4 
8.4    -2.874-10"4 
9.8    -3.001-10"4 
(STABLE) 

3.2    -7.459-10"4 
4.8    -1.096"10"3 
6.4    -9.78'10"4 
8.0    -1.121'10"3 

9.6    -1.534'10~3 

(STABLE  BUT 
INACCURATE) 

3.6    -8.327-10"4 
5.4    -1.982*10~3 
7.2    -2.335-10"3 
9.0    -3.801'10"3 
(UNSTABLE) 
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effect  of  the  Milne-predictor.   But  recalling  the  discussion 
of  the  published  results,  this  actual  real  negative  stability 
bound  is  exactly  the  same  as  that  obtained  by  Hamming  using 

two  iterations  with  truncation  error  modification.   This 

2 
coincidence  is  not  without  basis  since  the  P(EC)   mode  has 

exactly  the  same  truncation  error  modification  as  that  of 

the  Hamming  method,  making  the  two  experimental  procedures 

identical . 

6.  Second  Order  Adams  (P-C-VI) 

The  starting  value  chosen  was  h  =  0.1  with  an 
increment  of  0.1.   Table  7  shows  the  actual  data  obtained. 
For  h  =  0.1  up  to  h  =  0.3  the  method  is  stable.   But  for 
h  =  0.4  the  method  shows  actual  instability  as  evidenced  by 
the  uniformly  increasing  error  growth.   Thus  this  method 
is  stable  within  the  range  -0.4  <  hC  <  0.   Compared  to  the 
predicted  stability  of  the  corrector,  which  was  the  same  as 
Euler  CP-C-I),  with  range  of  stability  -2.0  <  hC  <  0,  it  is 
evident  that  the  stability  is  greatly  reduced.   Again  this 
could  be  traced  to  the  high  instability  of  the  predictor 
formula  used.   This  is  probably  the  reason  why  almost  all 
published  material  on  the  Adams -Moulton  method  always  starts 
with  the  third  order  and  higher  predictor.   Hull  and  Creemer 
[Ref.  8]  showed  that  this  method  has  the  largest  error  among 
the  Adams -type  methods. 

7.  Third  Order  Adams  (P-C-VII) 

The  choice  for  h  starts  at  h  =  0.8  and  is  then 
incremented  by  0.1.   Table  8  shows  the  selected  actual  data 
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TABLE  7 

ABSOLUTE  ERROR  (Z)  FOR  SECOND  ORDER  ADAMS  (P-C-VI) 

ODE:   y'  =  -y  with  y(0)=l 
Single  Precision 


h  =  0.1 

h  =  0.2 

x             ZA2 

x             ZA2 

1.0         3.017" 

10"4 

1.0         1.130-10"3 

2.0        2.862 

10'4 

2.0         1.347-10"3 

3.0         2.223 

10"4 

3.0         1.277-10"3 

4.0         1.772 

10'4 

4.0        1.196-10"3 

5.0         1.527 

io-4 

5.0        1.146-10"3 

6.0        1.408 

lO"4 

6.0        1.120-10"3 

7.0        1.354 

IO"4 

7.0        1.108-10"3 

8.0        1.330 

10"4 

-3 
8.0         1.103-10 

9.0        1.320 

■IO"4 

9.0        1.100-10"3 

10.0         1.315 

■IO"4 

10.0         1.099-10"3 

(STABLE) 

(STABLE) 

h  =  0.3 

h  =  0.4 

x             ZA2 

x             EA2 

1.2        2.607-10"3 

1.2        4.164-10"3 

2.4         3.589-10"3 

2.4        7.290-10"3 

3.6         3.814-10"3 

3.6        8.663-10"3 

4.8        3.861-10"3 

4.8        9.207-10"3 

6.0         3.865-10"3 

6.0        9.411-10"3 

7.2         3.869-10"3 

7.2         9.484-10"3 

8.4         3.868-10"3 

8.4        9.510-10"3 

9.6        3.868-10"3 

9.6        9.519-10"3 

(STABLE  BUT  INACCURATE) 

(UNSTABLE) 
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TABLE  8 

ABSOLUTE  ERROR  (E)  FOR  THIRD  ORDER  ADAMS  (P-C-VII) 

ODE:   y'=-y  with  y(0)=l 
Single  Precision 


h   =    0.8 

h    =    0.9 

x                                 EA3 

x                                 ZA3 

3.2                     1.402-10"2 
4.8                     7.751-10"3 
6.4                    1.968-10"3 
8.0                    4.614-10"4 
9.6                    4.608-10"4 
(STABLE) 

3.6                     1.600-10"2 
5.4                    8.272* 10"3 
7.2                     1.832-10"3 
9.0                    7.436-10"4 
(STABLE) 

h   =    1.2 

h    =    1.3 

x                                 ZA3 

x                                 EA3 

3.6                     7.06-10"5 
4.8                     8.773-10"3 
6.0                     5.019-10"3 
7.2                     5.216-10"3 
8.4                   -2.976-10"3 
9.6                     3.050*10"3 
(STABLE    BUT    INACCURATE) 

3.9                  -4.736-10"3 
5.2                     1.426'10"3 
6.5                     1.569-10"3 
7.8                    1.869'10"3 
(UNSTABLE) 
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obtained.   For  h  =  0.8  to  h  =  0.9  absolute  stability  occurs. 
For  h  =  1.0  to  h  =  1.2  oscillation  occurred  and  the  solution 
is  inaccurate  though  stable.   For  h  =  1.3  the  error  growth 
is  increasing  showing  actual  instability.   Thus  the  method 
is  stable  within  the  range  -1.3  <  hC  <  0 ,  which  is  in 
conformity  with  the  predicted  stability  of  the  general  Adams - 
method,  -1.3  <  hC  <  0 .   This  also  justifies  why  most  Adams  - 
Moulton  methods  start  with  the  use  of  the  third  order  predic- 
tor. 

8.   Fourth  Order  Adams  (P-C-VIII) 

Starting  value  was  h  =  0.5  with  0.1  increment. 
Table  9  shows  the  actual  data  obtained.   From  h  =  0.5  to 
h  =  0.7  the  method  is  absolutely  stable.   For  h  =  0.9  oscil- 
lation continues  as  with  h  =  0.8.   This  behavior  was  observed 
up  to  h  =  1.0  but  the  numerical  solution  is  still  stable 
though  inaccurate.   For  h  =  1.1  actual  instability  is  evident 
from  the  tendency  to  increasing  oscillation  of  the  error. 
Thus  this  method  is  stable  for  -1.1  <  hC  <  0 . 

Summarizing  then,  the  real  negative  stability  bounds  of 

the  different  P-C  sets  considered  are  outlined  below  with 

C  =  f  (x,y)  and  h  =  step  size. 

Real  Negative 
P-C  Set  Stability  Limit 

P-C-I  -1.9  <  hC  <  0 

P-C-II  Unstable 

P-C-III  -1.8  <  hC  <  0 

P-C-IV  Unstable 

P-C-V  -0.9  <  hC  <  0 

P-C-VI  -0.4  <  hC  <  0 

P-C-VII  -1.3  <  hC  <  0 

P-C-VIII  -1.1  <  hC  <  0 
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TABLE  9 

ABSOLUTE  ERROR  (E)  FOR  FOURTH  ORDER  ADAMS  (P-C-VIII) 

ODE:   y'  =  -y  with  y(0)=l 
Single  Precision 


h  =  0.5 

h  =  0.7 

x             ZA4 

x             EA4 

2.0        5.67-10"5 

3.0        2.125-10"4 

4.0        1.488-10"4 

5.0        8.11-10"5 

6.0         3.95-10"5 

7.0        1.8-10  "5 

8.0        7.9-10"6 

9.0        3.4-10"6 

10.0        1.4-10"6 
(STABLE) 

2.8        4.094-10"4 
4.2        8.157-10"4 
5.6        3.522-10"4 
7.0        1.239-10"4 
8.4         3.52-10"5 
9.8         7.7-10"6 
(STABLE) 

h  =  0.9 

h  =  1.1 

x             ZA4 

x             £A5 

3.6         1.819-10"3 
5.4        2.111*10"3 
7.2        4.93-10"5 
9.0       -5.125-10"4 
(STABLE  BUT  INACCURATE) 

3.3  -3.885-10"3 

4.4  5.789*10"3 

5.5  1.214-10"3 

6.6  3.818-10"3 

7.7  6.321-10"3 

8.8  -4.534-10"3 

9.9  7.846-10"3 
(UNSTABLE) 
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The  significant  points  noted  from  these  results  are: 

i)  Within  the  limits  of  their  stability,  the  Hamming 
and  Fourth  Order  Adams  methods  produced  the  best  accuracy. 

ii)  It  is  quite  interesting  to  see  that  although  the 
Euler  and  Nystrom  methods  have  the  widest  range  of  real 
negative  stability,  their  results  are  not  quite  as  accurae 
as  the  results  of  the  Hamming  and  the  Fourth  Order  Adams, 
even  the  Third  Order  Adams  methods. 

D.   NUMERICAL  RESULTS  FOR  RELATIVE  STABILITY 

If  C  >  0,  the  solution  itself  and  generally  also  the 
error  are  increasing  exponentially.   In  this  situation 
relative  stability  is  the  important  consideration.   A  P-C 
method  is  relatively  stable  if  the  rate  of  change  of  the 
error  with  respect  to  a  finite  range  of  integration  points 
is  less  than  the  rate  of  change  of  the  true  solution  with 
respect  to  the  same  finite  range  of  integration.   By  this 
definition  it  could  be  seen  that  it  is  extremely  difficult 
to  establish  a  fixed  bound  for  relative  stability.   In 
order  to  have  a  working  knowledge  of  the  relative  stability 
of  the  P-C  sets  considered  the  ODE 

y1  =  y  with  y(0)  =  1 

where  C  =  f  (x,y)  >  0  was  used.   Each  P-C  set  was  run  with 
the  same  series  of  h  values  from  h  =  0.1  to  h  =  2.0,  with 
an  increment  of  0.1   The  range  of  integration  for  each 
step  size  h  was  x  =  0  to  x  =  10.   Table  10  shows  the  true 
solution  values  from  x  =  1  to  x  =  10.0.   Table  10-A  to 
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TABLE  10 

TRUE  SOLUTION  VALUES  FOR 

ODE:   y'=y  with  y(0)=l 
Single  Precision 


^exact 


1.0  2.718 

2.0  7.389 

3.0  20.085 

4.0  54.597 

5.0  148.409 

6.0  403.417 

7.0  1096.595 

8.0  2980.838 

9.0  8102.710 

10.0  22025.320 
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Table  10-H  show  the  selected  actual  data  obtained  for  the 
different  P-C  sets. 

Table  10-A  shows  that  the  Euler  method  is  relatively 
stable  for  the  particular  range   of  integration.   Comparing 
the  error  growth  with  the  solution  growth  of  Table  10,  it 
could  be  seen  that  the  rate  of  change  of  error  is  less  than 
the  rate  of  change  of  the  solution.   From  h  =  0.1  to  h  =  0.5 
the  method  is  quite  accurate.  However  from  h  =  1.0  to 
h  =  2.0  it  becomes  quite  inaccurate  though  still  relatively 
stable.   Table  10-B  shows  that  the  Milne  method  is  not  only 
relatively  stable  but  accurate  from  h  =  0 . 1  to  h  =  2 . 0 . 
From  Table  10-C  the  Nystrom  method  exhibits  relative  stabi- 
lity from  h  =  0.1  to  h  =  1.0  but  becomes  unstable  at  h  =  2.0. 
This  could  be  determined  by  simply  comparing  the  magnitude 
of  the  error  at  x  =  10.0  and  h  =  2.0  with  the  magnitude  of 
the  true  solution  at  the  same  point,  in  which  case  it  is 
obvious  that  the  error  is  larger  than  the  true  solution. 
This  is  another  way  of  looking  for  relative  instability 
since  it  is  clear  that  if  the  rate  of  change  of  the  error 
is  greater  than  the  rate  of  change  of  the  solution  with 
respect  to  n  (the  number  of  solution  points)  then  eventually 
the  error  will  be  greater  than  the  true  solution.   Table  10-D 
shows  that  the  Hermite  method  is  relatively  stable  and  also 
accurate.   From  Table  10-E  it  can  be  seen  that  Hamming's 
method  is  also  stable  and  accurate.   Table  10-F  shows  that 
the  second  order  Adams'  method  is  stable  from  h  =  0.1  to 
h  =  1.0,  but  the  accuracy  is  quite  diminished  at  h  =  1.0. 
At  h  =  2.0  it  is  unstable.   From  Table  10-G,  the  third 
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TABLE  10 -A 

ABSOLUTE  ERROR  (Z)  FOR  EULER  (P-C-I) 

ODE:   y'=y  with  y(0)=l 
Single  Precision 


X 

h 

h   =    0.1 

h    =    0.5 

h    =    1.0 

h   =    2.0 

1.0 

2.0 
3.0 
4.0 
5.0 
6.0 
7.0 
8.0 
9.0 
10.0 

-2.231-10"3 

-1.219-10"2 

-4.980-10"2 

-1.806-10"1 

-6.142-10"1 

-2.004 

-6.360 

-19.774 

-60.496 

-182.828 

-4.819-10"2 

-2.570-10"2 

-1.046 

-3.808 

-13.011 

-42.721 

-136.450 

-427.070 

-1316.117 

-4006.574 

-1.567-10"1 

-7.515-10"1 

-2.959 

-10.638 

-36.261 

-119.355 

-383.286 

-1208.460 

-3756.500 

-11546.150 

-1.610 

-16.401 

-155.571 

-1420.042 

-12622.530 
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TABLE  10 -B 

ABSOLUTE  ERROR  (E)  FOR  MILNE  (P-C-II) 

ODE:   y'=y  with  y(0)=l 
Single  Precision 


X 

\ 

h  =  0.1 

h  =  0.5 

h  =  1.0 

h  =  2.0 

1.0 
2.0 
3.0 
4.0 
5.0 
6.0 
7.0 
8.0 
9.0 
10.0 

4.8-10"6 
-2.38-10"5 
-1.678-10"4 
-7.477-10"4 
-2.777-10"3 
-9.765-10"3 
-3.222- 10"2 
-1.040-10"1 
-3.281-10"1 
-1.003 

1.409*10"3 
-1.831*10"4 
-1.309-10"2 
-7.122-10"2 
-2.917-10"1 
-1.059 
-3.608 
-11.781 
-37.402 

3.070" 10"1 

8. 342-10"1 

1.794 

4.026 

8.289 
15.589 
23.242 

362.921 
3472.707 
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TABLE  10 -C 

ABSOLUTE  ERROR  (E)  FOR  NYSTROM  (P-C-III) 

ODE:   y'=y  with  y(0)=l 
Single  Precision 


X 

EN 

h  =  0.1 

h  =  0.5 

h  =  1.0 

h  =  2.0 

1.0 
2.0 
3.0 
4.0 
5.0 
6.0 
7.0 
8.0 
9.0 
10.0 

-2.034-10"3 

-1.173-10"2 

-4. 875  * 10-2 

-1.783-10"1 

-6.096-10"1 

-1.996 

-6.350 

-19.778 

-60.589 

-183.289 

-2.756-10"2 

-2.231-10"1 

-1.016 

-3.901 

-13.757 

-46.134 

-149.630 

-473.892 

-1474.339 

-4523.753 

-5.244-10"1 

-2.758 

-11.304 

-41.718 

-145. 106 

-485.911 

-1584.737 

-5069.113 

-15975.800 

-16.401 
-215.571 
-2384.042 
-24472.530 
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TABLE  10 -D 

ABSOLUTE  ERROR  (Z)  FOR  HERMITE  (P-C-IV) 

ODE:   y'=y  with  y(0)=l 
Single  Precision 


X 

h  =  0.1 

h  =  0.5 

h  =  1.0 

h  =  2.0 

1.0 

2.0 
3.0 
4.0 
5.0 
6.0 
7.0 
8.0 
9.0 
10.0 

3.8-10"6 
-2. 77-10"5 
-1.831-10"4 

-7.629-10"4 
-2.853-10"3 
-9.765-10"3 
-3.247-10"2 
-1.044-10"1 
-3.281-10"1 
-1.003 

-3.948-10"4 

-2.104-10"3 

-9.292-10"3 

-3.546-10"2 

-1.244-10"1 

-4.147-10"1 

-1.335  ' 

-4.194 

-12.941 

-39.367 

-1.217-10"2 

-2.285-10"2 

-8.180-10"2 

-2.519-10"1 

-7.766-10"1 

-2.355 

-7.-67 

-21.023 

-62.078 

7.092-10"1 
21.334 
275.2517 
2860.671 
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TABLE  10 -E 

ABSOLUTE  ERROR  (E)  FOR  HAMMING  (P-C-V) 

ODE:   y'=y  with  y(0)=l 
Single  Precision 


X 

EHA 

h  =  0.1 

h  =  0.5 

h  =  1.0 

h  =  2.0 

1.0 
2.0 
3.0 
4.0 
5.0 
6.0 
7.0 
8.0 
9.0 
10.0 

2.19-10"5 

6. 77-10"5 

1.831-10"5 

6.714-10"4 

1.831-10"3 

5.371-10"3 

1.879- 10~2 

4.687*10"2 

1.406-10"1 

4.296-10"1 

2.980-10"3 
2.792-10"3 

-7.202-10"3 

-6.004-10"2 

-2.777- 10"1 

-1.039 

-3.639 

-12.097 

-38.882 

4.332-10"1 
8.498-10"1 
1.043 

-5.24-10"1 

-10.532 

-53.558 

-213.566 

380.836 
3040.449 
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TABLE  10 -F 

ABSOLUTE  ERROR  (E)  FOR  SECOND  ORDER  ADAMS  (P-C-VI) 

ODE:   y'=y  with  y(0)=l 
Single  Precision 


X 

h 

h   =    0.1 

h    =    0.5 

h   =    1.0 

h   =    2.0 

1.0 
2.0 
3.0 
4.0 
5.0 
6.0 
7.0 
8.0 
9.0 
10.0 

-2.100-10"3 

-1.242-10"2 

-5.241-10"2 

-1.934-10"1 

-6.649-10"1 

-2.186 

-6.974 

-21.766 

-66.800 

-202.378 

-2.735-10"2 

-2.358-10"1 

-1.127 

-4.453 

-16.013 

-54.444 

-178.434 

-569.834 

-1785.156 

-5510.332 

-5.041-10"1 

-2.792 

-11.936 

-45.345 

-161.08 

-548.354 

-1812.745 

-5866.652 

-18684.180 

-14.40.1 

-201.571 

-2324.042 

-24506.530 
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TABLE  10 -G 

ABSOLUTE  ERROR  (E)  FOR  THIRD  ORDER  ADAMS  (P-C-VII) 

ODE:   y'=y  with  y(0)=l 
Single  Precision 


X 

EA3 

h  =  0.1 

h  =  0.5 

h  -  1.0 

h  =  2.0 

1.0 

2.0 
3.0 
4.0 
5.0 
6.0 
7.0 
8.0 
9.0 
10.0 

2.155-10"4 

1.318-10"3 

5.569-10"3 

2.062-10"2 

7.075-10"2 

2.324-10"1 

7.412-10"1 

2.307 

7.089 

21.468 

2.076-10"2 

1.454-10"1 

6.352-10"1 

2.376 

8.216 

27.090 

86.509 

269.992 

828.183 

2.131-10"2 

5.538-10"1 

2.827 

11.147 

39.645 

132.954 

429.273 

1349. 750 

14.938 

217.016 

2409.128 
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TABLE  10- H 

ABSOLUTE  ERROR  (Z)  FOR  FOURTH  ORDER  ADAMS  (P-C-VIII) 

ODE:   y'=y  with  y(0)=l 
Single  Precision 


X 

EA4 

h   =    0.1 

h   =    0.5 

h   =    1.0 

h   =    2.0 

1.0 
2.0 
3.0 
4.0 
5.0 
6.0 
7.0 
8.0 
9.0 
10.0 

8.6-10"6 
-2.77-10"5 
-1.984-10"4 
-8.087-10"4 
-3.234-10"3 
-1.074*10"2 
-3.540-10"2 
-1.176-10"1 
-3.476-10"1 
-1.085 

S.53-10"5 
-1.995-10"2 
-1.089-10"1 
-4.450-10"1 
-1.614 
-5.490 
-17.926 
-56.878 
-176.812 

9.992-10"2 
-7.852-10"1 
-5.094 
-21.956 
-81.859 
-2  83.2  30 
-936.164 

262.678 
1748.47 

■     - 
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order  Adams'  method  showed  accuracy  and  stability  from 
h  =  0.1  to  h  =  2.0.   Finally,  Table  10-H  shows  that  the 
Fourth  order  Adams'  method  is  quite  stable  and  accurate 
from  h  =  0.1  to  h  =  2.0. 

The  most  interesting  points  observed  from  these  actual 
results  are: 

i)   Though  the  Milne  and  Hermite  methods  have  no  real 
negative  stability  limits,  they  are  quite  accurate  for 
C  >  0.   In  fact  the  Hermite  method  produced  the  least  start- 
ing error,  and  the  Milne's  method  provided  the  second  least 
starting  error,  and  both  maintained  accuracy  up  to  h  =  2.0. 

ii)   In  contrast,  the  Euler  and  Nystrom  methods  both 
have  the  widest  range  of  real  negative  stability  limits  but 
showed  inaccurate  results  starting  at  h  =  1.0.   In  fact  the 
Nystrom  method  is  unstable  at  h  =  2.0. 

iii)   The  third  and  fourth  order  Adams'  and  Hamming  methods 
showed  wide  range  of  relative  stability  and  provided  accurate 
results.   The  second  order  Adams'  method  again  showed  inac- 
curate results  at  h  =  1.0  and  relative  instability  at  h  =  2.0. 

iv)   From  h  =  0.1  to  h  =  1.0,  the  methods  in  terms  of 
accuracy,  rank  as  follows:   Milne  ranks  first,  Hermite  second, 
then  Hamming,  Fourth  Order  Adams,  and  Third  Order  Adams  in 
that  order. 
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IX.   NUMERICAL  EXPERIMENTS 

In  order  to  test  the  performance  of  the  different  P-C 
sets  considered,  a  collection  of  test  ODEs  with  many  unusual 
and  interesting  features  (i.e.,  singularities,  discontinui- 
ties, infinite  derivatives,  oscillating  derivatives,  etc.) 
were  selected.   Several  of  such  ODEs  were  presented  by  Hull 
and  Creemer  [Ref .  8]  and  Lapidus  and  Seinfeld  [Ref .  2] . 
Table  11  lists  the  test  ODEs  considered.   To  simplify  nota- 
tion, the  ODE  number  will  be  used  to  specify  the  differential 
equation.   For  example  ODE  I  is  equivalent  to  the  ODE 
y'  =  -y  +  10sin3x  with  y(0)  =  -3,  whose  analytic  solution  is 
y(x)  =  sin3x  -  3cos3x. 

From  the  previous  section  on  stability  bounds  of  P-C 
sets  it  was  observed  that  the  different  methods  showed  good 
accuracy  up  to  h  =  0.5  on  both  experimental  ODEs,  y'  =  -y 
and  y'  =  y,  considered.   It  was  further  observed  from  pre- 
vious numerical  results  that  all  the  P-C  sets  exhibited 
good  stability  behavior  within  the  range  of  h  up  to  h  =  0.5 
except  for  the  Milne  and  Hermite  methods,  in  the  case  of 
real  negative  stability.   But  these  shortcomings  of  the  Milne 
and  Hermite  methods  were  compensated  by  the  fact  that  they 
produced  good  accuracy  and  wide  range  of  relative  stability. 
This  analysis  is  needed  in  the  sense  that  the  test  ODEs  con- 
sidered exhibit  both  cases  of  f  (x,y)  <  0  and  f  (x,y)  >  0. 
By  using  values  of  h  =  a   ,2   ,2   ,2   ,2   ,2   ,2 
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where  the  largest  value  of  h  =  0.5,  each  method  will  be 
expected  to  perform  quite  well.   Thus  a  comparative  analysis 
of  their  accuracy  and  computing  time  would  be  meaningful. 
In  dealing  with  accuracy,  the  effects  of  the  propagation 
of  errors  would  be  strongly  felt,  as  has  already  been  shown 
in  the  previous  analysis  of  error  propagation.   While 
computing  time  takes  into  its  fold  the  number  and  complexity 
of  function  evaluations. 

All  the  eight  P-C  sets  were  run  on  each  test  ODE  using 
the  series  of  h  values  as  previously  mentioned.   The  Fortran 
programs  used  were  the  same  as  those  used  for  stability 
analysis,  but  the  precision  used  was  Double  Precision  (14 
digit  precision  for  IBM  360/67  computer)  to  minimize  the 
effect  of  rounding  error  as  h  decreases.   Each  test  ODE  will 
be  studied  individually,  with  the  objective  of  providing 
useful  and  important  recommendations  on  which  P-C  set 
performs  best  for  the  particular  class  of  problems. 

A.   ODE  I 

Table  12  shows  the  exact  solution  from  x  =  1.0  to 
x  =  10.0.   From  Table  12-A  the  results  show  that  the  Euler 
method  provides  the  best  accuracy  with  7.65  sees  in  computing 
time  with  the  Nystrom's  ranking  second  in  accuracy  though 
with  less  computing  time,  7.46  sees.   The  Milne  and  the 
Hermite  methods  are  not  considered  since  they  both  showed 
instability  in  their  numerical  solution  values.   This  could 
easily  be  seen  by  looking  at  Table  12-B  which  shows  the 
behavior  of  the  errors  for  these  methods.   The  errors  for 
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TABLE  12 

TRUE  SOLUTION  VALUES  FOR  ODE  I 
Double  Precision 


^exact 


1.0  3.111097 

2.0  -3.159926 

3.0  3.145509 

4.0  -3.068134 

5.0  2.929351 

6.0  -2.731937 

7.0  2.479843 

8.0  -2.178115 

9.0  1.832792 

10.0  -1.450785 
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TABLE  12-A 

ODE  I  AT  x  =  10.0 

ABSOLUTE  ERROR  FOR  P-C-I,  P-C-II,  P-C-III,  and  P-C-III 

Double  Precision 


v 


1/h 

h 

M 

h 

^HE 

128 

64 

32 

16 

8 

4 

2 

CPTS* 

-9.83479-10"5 
-3.93360-10"4 
-1.57286-10"3 
-6.28126-10"3 
-2.49436-10"2 
-9.64382-10"2 
-3.34570-10"1 
7.65 

4.51483-10"9 

8.21731-10"8 

-1.17807-10"7 

-2.51544-10"5 

-1.35324-10"3 

-3.35324-10"2 

1.527752 

5.59 

-9.83522-10"5 
-3.93415'10"4 
-1 . 57372 ' 10"3 
-6.92664-10"3 
-2.52237-10"2 
-1.00786-10"1 
-3.38212-10"1 
7.46 

-8.53288-10"3 
-3.40492-10"2 
-1.35592,10~1 
-5.38364-10"1 
-2.134035 
-6.039341 
-37.089294 
8.37 

CPTS  means  computing  time  in  seconds. 
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TABLE  12 -B 

ABSOLUTE  ERROR  FOR  P-C-I,  P-C-II,  P-C-III,  and  P-C-IV 

Double  Precision 
ODE  I ;  h  =  0.5 


'M 


'N 


'HE 


1.0 

2.0 
3.0 
4.0 
5.0 
6.0 
7.0 
8.0 
9.0 
10.0 


7.171-10 


4.640-10 


-1 


5.756'10 


-5.392-10 


5.418-10 


5.154-10 


4.854*10 


-4.431*10 


-1 


-1 


-1 


3.929-10"1 


3.345 


2.081*10 


9.447-10 


2.509-10"1 


2.076-10 

4.266-10 

4.737-10 

7.910*10 

1.011 

1.527 


4.159-10 


-6.649-10 


5.817-10 


6.091-10 


-1 


-1 


5.806*10 


-5.560*10 


5.148-10 


-4.652*10 


4.056*10 


3.382 -10 


1.584 

1.645 

2.464 

3.386 

5.218 

7.537 

11.400 

11.747 

25.066 

37.089 
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the  Milne  and  Hermite  methods  increased  in  one  direction 
eventually  growing  to  magnitudes  greater  than  the  true 
solution  at  x  =  10.0  (i.e.,  for  Milne's  error,  1.527  >  1.450 
and  for  Hermite's  error  -37.089  >  1.450),  thus  these  methods 
are  not  suited  for  solving  ODE  I.   The  Euler  and  Nystrom 
methods  both  started  at  accuracy  approximately  equal  to 
fifth  decimal  places  and  ended  up  at  h  =  0 . 5  with  accuracy 
up  to  the  first  decimal  place.   But  though  their  accuracy 
is  quite  restrained  as  h  increases,  nevertheless  the  behavior 
of  the  error  tends  to  decrease  as  the  range  of  integration 
is  increased.   From  this  group  then  the  Euler  method  is 
the  best  choice  for  ODE  I. 

In  Table  12-C  the  results  for  the  other  four  methods  are 
listed.   All  these  methods  showed  stability  for  all  values 
of  h  used,  but  Hamming's  produced  the  best  accuracy  and 

the  best  computing  time.   Hamming's  method  started  at  10 

_  2 
accuracy  at  1/h  =  128  and  ended  with  10    at  h  =  0.5.   The 

Fourth  Order  Adams'  comes  next  in  both  accuracy  and  computing 

time,  then  the  Third  Order  Adams',  and  the  Second  Order 

Adams'  came  in  that  order. 

Comparing  these  last  four  methods  with  the  Euler  method, 

only  Adams'  second  order  method  is  inferior.   Therefore,  for 

solving  ODEs  that  belong  to  the  class  of  ODE  I  the  Hamming 

method  is  highly  recommended  in  terms  of  accuracy  and  least 

cost  in  computer  time.   The  Fourth  Order  Adams'  is  the  next 

choice  in  the  absence  of  the  Hamming  method. 
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TABLE  12 -C 

ODE  I  at  x  =  10.0 

ABSOLUTE  ERROR  (£)  FOR  P-C-V,  P-C-VI,  P-C-VII,  and  P-C-VIII 

Double  Precision 


1/h 

^Ha 

ZA2 

ZA3 

ZA4 

128 

64 

32 

16 

8 

4 

2 

CPTS 

3.0381-10"10 
1.12626-10"8 
4.68494-10"7 
1.01544'10"5 
3.69762-10"4 
1.29180-10"2 
6.39906-10"2 
5.76 

-9.86551-10"5 
-3.95790-10"4 
-1.59200-10"3 
-6.42959-10"3 
-2.60335-10"2 
-1.03170-10"1 
-3.586413 
7.98 

-9.60790*10"7 
-7.81167-10"6 
-6.45309-10"5 
-5.49889-10"4 
-4.97320-10"3 
-5.00046-10"2 
-4.74203-10"1 
8.55 

1.86070-10"8 
3.18596-10"7 
4.04497-10"6 
5.91417-10"5 
7.48403-10"4 
4.56055'10"3 
-1.90291-10"1 
5.99 
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B.   ODE  II 

The  true  solution  values  are  given  in  Table  13  from 
x  =  1.0  to  x  =  10.0.   From  the  data  presented  in  Table  13-A 
it  is  clear  that  the  most  important  consideration  is  the 
step  size  h.   For  h  =  0.5  the  four  methods  identically 
yield  errors  greater  than  the  true  solution  at  the  same 
point.   In  solving  problems  of  this  type  then  h  must  be 
chosen  to  be  smaller  than  0.25  for  the  numerical  solution 
to  be  valid  using  the  Euler,  Milne,  and  Nystrom  methods. 
The  Hermite  method  obviously  is  not  worth  considering, 
because  of  its  inaccurate  results  for  all  values  of  h. 
For  h  <0.25  the  Milne  method  gives  the  best  numerical 
solution  and  the  least  computing  time,  followed  by  the  Euler 
then  the  Nystrom  methods. 

Table  13-B  shows  that  to  have  a  vlid  numerical  solution 
for  any  of  the  other  four  methods  h  must  be  less  than  0.25. 
For  values  of  h  <  0.25  ,  the  Hamming  method  yields  the 
greatest  accuracy,  followed  by  the  Fourth,  Second,  and 
Third  Order  Adams'.   The  Milne  method  is  a  little  better 
than  the  Second  Order  Adams'.   Thus  for  accuracy  the  Hamming 
method  ranks  first,  then  the  Fourth  Order  Adams,  Milne, 
Euler,  Second  Order  Adams,  Third  Order  Adams,  and  lastly 
the  Nystrom  methods. 

For  problems  in  the  class  of  ODE  it  is  recommended  that 
h  must  be  chosen  less  than  0.25  or  smaller  if  valid  numerical 
solution  and  accuracy  are  desired.   Then  use  Hamming  as  the 
numerical  method  to  solve  the  ODE.   Again  the  Fourth  Order 
Adams  method  is  the  second  best  choice.   The  choice  of  h  can 
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TABLE  13 

TRUE  SOLUTION  VALUES  FOR  ODE  II 
Double  Precision 


y 

7  exact 


1.0  -1.381773 

2.0  -4.931505-10"1 

3.0  8.488724-10"1 

4.0  1.410446 


5.0  6.752620-10"1 


6.0  -6.807547-10"1 


7.0  -1.410888 

8.0  -8.438582-10"1 

9.0  4.990117-10"1 

10.0  1.383092 
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TABLE  13-A 

ODE  II  at  x  =  10.0 

ABSOLUTE  ERROR  FOR  P-C-I,  P-C-II,  P-C-III,  and  P-C-IV 

Double  Precision 


1/h 


*M 


JN 


'HE 


128 
64 
32 
16 


2.596512-10 


-3 


2.50786*10 


-9 


2.518638-10 


2.726949*10 


-3 


4.095807-10 


1.634916*10 

9.921730-10" 
-4 


CPTS 


6.213302-10 


9.160642 


131.604044 


6.38 


-1 


5.69316710 


2.178410-10 


1.362692-10 


5.482274 


4.69 


-3 


-1 


8.444203*10 
7.950519-10 
5.455366-10 
4.234564-10 
3.214692 
23.407199 
157.762746 
5.34 


-4.446340-10 
-1.764591 
-6.947361 
-26.911798 
-100.834152 
-353.009679 
-1081.991963 
8.26 
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TABLE  13-B 

ODE  at  x  =  10.0 

ABSOLUTE  ERROR  FOR  P-C-V,  P-C-VI,  P-C-VII,  and  P-C-VIII 

Double  Precision 


1/h 


'HA 


JA2 


'A3 


*A4 


128 

64 

32 

16 

8 

4 


1.20446*10 


7.48373- 10 


-9 


-9 


3.760418*10 


-4 


-3.192068-10 


5.531148-10 


-2.331998-10 


7.353901-10 


3.686296-10 


-3 


1.532622-10 


7.673571-10 


-1 


-1 


1.657253-10 


6.292332 


CPTS 


5.26 


-6. 706531-10 


44.871662 


7.79 


4.448046-10 


3.610822-10 


2.961989*10 


2.455737-10 


2.016912 


15.053093 


80.187459 


8.36 


2.334506-10 


-6.282905-10 


-1.347495-10 


-6.597818-10 


-1.600424-10 


-2 


2.268373-10 


1.401846 


4.96 
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be  formally  stated  as  |hC|  <  0.25,  since  the  value  of  h  to 
be  used  is  directly  dependent  upon  C  =  f  (x,y). 

C.  ODE  III 

Table  14  shows  the  exact  solution  values  for  ODE  III  from 
x  =  1.0  to  x  =  10.0.   From  the  data  of  Table  14-A  it  is  clear 
that  all  four  methods  performed  quite  well  to  a  good  degree 
of  accuracy  for  all  balues  of  h.   But  for  best  accuracy  and 
least  computing  time  the  Milne  Method  stands  out  followed  by 
the  Euler,  Hermite,  and  Nystrom  methods  in  that  order. 

From  Table  14-B,  again  it  is  observed  that  any  of  the  last 
four  methods  yield  a  valid  numerical  solution,  to  a  certain 
degree  of  accuracy.   For  an  excellent  accuracy,  however,  the 
Hamming  method  is  the  best  choice.   Its  range  of  accuracy  is 
from  10    up  to  10    for  values  of  h  ranging  from  1/h  =  128 
to  1/h  =  2.   The  Milne  method  compared  to  these  last  four 
methods  ranks  second  only  to  the  Hamming,  though  its  com- 
puting time  is  a  little  smaller  than  Hamming's.   Thus  for 
excellent  accuracy,  the  order  of  choice  is  an  follows: 
Hamming,  Milne,  Fourth  Order  Adams,  Third  Order  Adams,  Her- 
mite, Euler,  Second  Order  Adams,  and  lastly  the  Nystrom 
method  in  solving  ODEs  belonging  to  the  class  of  ODE  III. 

D.  ODE  IV 

True  solution  values  are  listed  in  Table  15  from  x  =  1.0 
to  x  =  10.0.   Actual  data  obtained  in  Table  15-A  showed  that 
all  four  P-C  sets  are  good  numerical  methods  for  ODE  IV.   But 
if  a  choice  is  to  be  made  the  Milne's  accuracy  is  far  greater 
than  the  other  three  methods  with  computing  time  only  0.16 
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TABLE  14 

TRUE  SOLUTION  VALUES  FOR  ODE  III 
Double  Precision 


^exact 


1.0  8.414709-10"1 

2.0  9.092974-10"1 

3.0  1.411200-10"1 

4.0  -7.568024-10"1 

5.0  -9.589242-10"1 

6.0  -2.794154-10-1 

7.0  6.569865-10"1 

8.0  9.893582-10"1 

9.0  4.121184-10"1 

10.0  -5.440211-10"1 
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TABLE  14-A 

ABSOLUTE  ERROR  FOR  P-C-I,  P-C-II,  P-C-III,  and  P-C-IV 

ODE  III  AT  x=10.0 
Double  Precision 


1/h 


'M 


'N 


'HE 


128 

64 

32 

16 

8 

4 

2 

CPTS 


-1.705737*10 


-6.845247-10 


-2.703738-10 


-5 


-1.089843-10 


4.365210-10 


-4 


1.755574-10 


7.174787-10 


7.81 


-3 


1.099-10 


9.286-10 


11 


11 


1.785639-10 


1.341400-10 


1.26041-10 


-9 


1.58770-10 


1.596440-10 


7.814645-10 


3.531438-10 


1.736454-10 


-6 


-5 


-4 


1.288789-10 


1.031052-10 


8.247845-10 


-5 


-5 


-6 


4. 317133* 10 


1.071599-10 


5.59 


-3 


-9.528314-10 


5.856889-10 


3.910558-10 


7.12 


6.594418-10 


5.260049-10 


4.154229-10 


7.92 


-2 
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TABLE  14 -B 

ODE  AT  x  =  10.0 

ABSOLUTE  ERROR  FOR  P-C-V,  P-C-VI,  P-C-VII,  and  P-C-VIII 

Double  Precision 


1/h 


'HA 


*A2 


'A3 


'A4 


128 

64 

32 

16 

8 

4 

2 

CPTS 


2.04*10 


12 


6.396-10 


-11 


2.01370*10 


-9 


6.236927'IQ 


1.970409-10 


6.368642*10 


1.569286-10 


6.23 


5.023475-10 


-7.798238-10 


3.512798'IQ 


1. 721024-10 


1.900550-10 


6.720151*10 


-8 


-7 


3.683-10 


11 


6.5237-10 


10 


-5.373282*10 


4.287630*10 


9.408691-10 


5.767408-10 


-4 


-3 


3.391988-10 


-4 


3.848313- 10 


10.0 


2.588546-10 


-1.686135-10 


1.246830-10 


2.649671-10 


6.946177-10 


1.805833-10 


10.0 


4.360227-10 


8.23 


-3 
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TABLE  15 

TRUE  SOLUTION  VALUES  FOR  ODE  IV 
Double  Precision 


y 

}  exact 


1.0  3.894003 

2.0  3.639183 

3.0  3.306565 

4.0  2.943035 

5.0  2.578543 

6.0  2.231301 

7.0  1.911513 

8.0  1.624023 

9.0  1.370189 

10.0  1.149189 
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TABLE  15-A 

ODE  IV  at  x  =  10.0 

ABSOLUTE  ERROR  FOR  P-C-I,  P-C-II,  P-C-III,  and  P-C-IV 

Double  Precision 


1/h 


JM 


'N 


'HE 


128 

64 

32 

16 

8 

4 

2 

CPTS 


-7.344900-10 


-7 


-1.839350-10 


3.4*10 


6.07-10 


4.598339*10 


-12 


1.840047-10 


7.335362-10 


-6 


1.0603-10 


10 


7.339871-10 


-6 


-2.926381-10 


-1.176812-10 


-4 


1.97060-10 


-4.063982-10 


-8 


4.709013-10 


-1.885932- 10 


-6.024899-10 


4.39 


-5.791009*10 


3.83 


-6 


2.872267*10 


-1.119506-10 


4.260573-10 


1.541727-10 


3.67 


4.575814-10 


1.826764-10 


7.278726*10 


2.889080-10 


1.128099-10 


-6 


-5 


-5 


-4 


-3 


4.418111-10 


1.668370-10 


6.74 
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sec  longer  than  the  next  best  choice  which  is  Nystrom.   After 
Nystrom,  the  Euler  method  is  preferred,  then  the  Hermite. 

Table  15-B  data  again  demonstrates  that  any  of  the  last 
four  P-C  sets  can  be  used  for  solving  ODE  IV.   Milne  method 
compared  with  these  methods  competes  with  the  Hamming  in 
accuracy  and  has  an  edge  in  computer  time  (3.83  sees  <  4.45 
sees) .   But  this  fraction  of  seconds  difference  is  overcome 

by  the  demonstrated  better  accuracy  of  the  Hamming  in  every 

-  7         -1 
step  of  the  solution  from  h  =  2    to  h  =  2   .   Thus,  if 

preference  is  to  be  made,  Hamming  is  the  most  likely  choice 

as  the  best  method  for  solving  ODE  IV  class  of  problems. 

The  Milne,  Fourth  Order  Adams,  Third  Order  Adams,  Nystrom, 

Second  Order  Adams,  Euler,  and  Hermite  are  the  order  of 

choices  following  the  Hamming  method. 

E.   ODE  V 

Table  16  shows  the  true  solution  values  for  the  range  of 
integration  considered.   From  the  data  of  Tables  16-A  and 
16-B,  it  is  obvious  that  any  method  can  be  used  to  solve 
ODE  V  and  obtained  accuracy  up  to  minimum  of  10    and 
maximum  of  10    .   Thus  in  comparing  these  methods,  accuracy 
criteria  must  not  be  the  greatest  concern.   By  studying 
closely  the  behavior  of  the  errors  as  h  increases,  it  was 

noted  that  the  Milne,  Hamming,  and  Fourth  Order  Adams 

-7  -5 

exhibited  decreasing  errors  from  h  =  2    to  h  =  2    for 

Milne,  and  from  h  =  2    to  h  =  2    for  both  Hamming  and 

Fourth  Order  Adams.   Analyzing  further  the  Hamming  and 

Fourth  Order  Adams  methods  it  could  be  seen  that  the  Fourth 
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TABLE  15 -B 

ODE  IV  at  x  =  10.0 

ABSOLUTE  ERROR  FOR  P-C-V,  P-C-VI,  P-C-VII,  and  P-C-VIII 

Double  Precision 


1/h 


'HA 


'A2 


'A3 


'A4 


128 

64 

32 

16 

8 

4 

2 

CPTS 


3.0-10 


2.4-10 


13 


13 


2.37-10 


-12 


7.849-10 


11 


4.585327-10 


1.829584-10 


7.282200-10 


-6 


-6 


2.78944-10 


-9 


1.186192-10 


2.202890-10 


4.45 


-6 


2.884380*10 


-1.131593-10 


4.359162*10 


-1.625329*10 


6.46 


1.28499-10 


1.025906-10 


8.174763-10 


6. 488701*10 


1.61*10 


12 


2.690-10 


11 


4.3747-10 


-10 


5.111696*10 


7.20867*10 


-1.225493-10 


3.970077*10 


-5 


1.778276*10 


-6 


3.006535-10 


7.02 


2.201592*10 


4.02 
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TABLE  16 

TRUE  SOLUTION  VALUES  FOR  ODE  V 
Double  Precision 


y 

J  exact 


1.0  1.020204 

2.0  1.040832 

3.0  1.061913 

4.0  1.083472 

5.0  1.105541 

6.0  1.128152 

7.0  1.151338 

8.0  1.175139 

9.0  1.199593 

10.0  1.224744 
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TABLE  16 -A 

ODE  V  at  x  =  10.0 

ABSOLUTE  ERROR  FOR  P-C-I,  P-C-II,  P-C-III,  and  P-C-IV 

Double  Precision 


1/h 


*M 


*N 


'HE 


128 

64 

32 

16 

8 

4 

2 

CPTS 


2.26254*10 


-9 


9.04072-10 


3.609166-10 


3.436058-10 


7.0-10 


4.0-10 


14 


14 


-1.94002-10 


-9 


-4.8-10 


13 


7.747767*10 


-4.45-10 


-12 


0* 


2.2-10 


-3.088771-10 


3.583*10 


11 


13 


1.227293-10 


-7 


2.8689-10 


-10 


4.960094-10 


1.978142-10 


-7 


-6 


3.28-10 


12 


4.237-10 


11 


7.866820-10 


5.63 


-4.0408-10 


5.58 


10 


4.844143-10 


1.886736*10 


7.644850-10 


5.16 


-2.29641*10 


-1.839847-10 


1.478586-10 


5.35 


-7 


*  - 14 

0  means  error  less  than  10 
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TABLE  16-B 

ODE  V  at  x  =  10.0 

ABSOLUTE  ERROR  FOR  P-C-V,  P-C-VI,  P-C-VII,  and  P-C-VIII 

Double  Precision 


1/h 


'HA 


'A2 


'A3 


'A4 


128 
64 
32 
16 


4.0-10 


13 


1.9-10 


1.2-10 


5.0*10 


1.4*10 


-13 


-13 


-14 


-13 


4 
2 
CPTS 


-4.83*10 


12 


1.1744-10 


6.31 


10 


9.29021*10 


-9 


6.0-10 


14 


1.6-10 


-13 


1.179555-10 


2.46-10 


12 


1.1-10 


3.331471*10 


-8 


2.018-10 


11 


7.0-10 


13 


14 


1.249903*10 


4.926236-10 


1.950214*10 


-1.6114*10 


10 


-3.0*10 


-14 


-1.27940-10 


-1.007531-10 


1.14*10 


12 


-8 


1.802*10 


11 


7.649028*10 


10.0 


-6 


-7.798880-10 


10.0 


2.7389-10 


7.78 


-10 
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-4 
Order  Adams  had  its  optimum  value  of-h  at  h  <  2    while  that 

-3 
of  Hamming  achieved  its  optimum  h  value  at  h  <  2   .   Thus 

the  Hamming  method  has  a  higher  value  of  h  for  its  range  of 
excellent  accuracy.   This  property  gives  the  Hamming  method 
a  slight  edge  over  that  of  Milne,  and  the  Fourth  Order  Adams 
However  if  both  accuracy  and  computing  time  are  considered 
the  Milne  method  is  likely  to  be  the  choice.   Therefore  for 
the  class  of  ODE  V,  the  recommendation  for  the  best  method 
to  be  used  will  be  most  likely  dependent  upon  the  particular 
interest:   whether  accuracy  and  wider  range  of  h  values  are 
the  primary  concern,  or  whether  accuracy  and  least  computing 
time  is  the  criterion.   For  the  former  criteria  the  Hamming 
method  serves  the  best  purpose  while  for  the  latter  the 
Milne  method  offers  the  best  solution.   However,  as  was 
noted  earlier,  if  no  criterion  is  involved  but  the  interest 
is  just  to  solve  the  problem  the  easiest  way,  the  self- 
starting  Euler  method  is  recommended.   For  realistic  pur- 
poses however  the  following  order  of  choice  is  highly 
recommended:   Hamming,  Milne,  Fourth  Order  Adams,  Third 
Order  Adams,  Hermite,  Nystrom,  Euler,  and  lastly  Second 
Order  Adams . 

F.   ODE  VI 

True  solution  values  are  listed  in  Table  17  from  x  =  1.0 
to  x  =  10.0.   Results  from  Table  17-A  clearly  showed  that  h 
must  be  chosen  to  be  quite  small  to  have  a  valid  numerical 
solution.   For  the  series  of  h  values  chosen  none  of  the 
methods  exhibited  valid  numerical  solution  from  1/h  =  32  up 
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TABLE  17 

TRUE  SOLUTION  VALUES  FOR  ODE  VI 
Double  Precision 


^exact 


1.0  1.732050 

2.0  2.236067 

3.0  2.645751 

4.0  3.0 

5.0  3.316624 

6.0  3.605551 

7.0  3.872983 

8.0  4.123105 

9.0  4.358898 

10.0  4.582575 
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TABLE  17- A 

ODE  VI  at  x  =  10.0 

ABSOLUTE  ERROR  FOR  P-C-I,  P-C-II,  P-C-III,  and  P-C-IV 

Double  Precision 


1/h 

h 

EM 

ZN 

^HE 

128 

64 

32 

16 

8 

4 

2 

CPTS 

-50.844 
-107.099 
-214.685 
-436.640 
-898.232 
-1940.153 
-4695.585 
5.65 

-9.503-10"3 
7.761'10"2 
6.986 

-12.200 

-55.310 

-210.374 

-879.472 
4.40 

-47.611 
-99.681 
-198.306 
-373.796 
-660.199 
-1051.118 
-1475.333 
5.12 

-32.764 
-1.885 
18518.873 
382.107 
-3060.609 
-2279.396 
-4096.870 
7.11 
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to  h  =  2   .   The  Milne  method  provided  a  valid  numerical 
solution  starting  at  h  =  2    and  smaller,  but  the  other 
methods  still  are  unreliable. 

From  Table  17-B,  only  the  Hamming  and  Fourth  Order  Adams 
method  showed  a  valid  numerical  solution  for  h    2    and 

smaller,  all  others  need  much  smaller  h  than  the  lowest 

_  7 
value,  2    chosen.   Again  on  the  average  comparing  the  Milne, 

Hamming,  and  the  Fourth  Order  Adams,  for  use  in  solving 

ODE  VI  with  values  of  h    2   ,  the  Hamming  is  the  first 

choice,   Milne  second  and  the  Fourth  Order  Adams.   As  in 

the  case  of  ODE  II,  the  h  values  are  governed  by  the  maximum 

magnitude  of  C  =  f  (x,y).   If  C  is  large  then  h  must  be 

chosen  to  be  very  small  to  satisfy  the  bounds  for  stability 

as  had  been  established  before. 


G.   ODE  VII 

True  solution  values  are  shown  in  Table  18  from  x  =  1.0 
to  x  =  10.0.   Table  18-A  showed  that  the  Euler  and  Nystrom 
provide  the  accurate  numerical  solutions  desired  while  the 
Milne  and  Hermite  are  unreliable.   By  analyzing  the  behavior 
of  the  roots  of  the  Milne  and  Hermite  methods  shown  in 
Table  18-B,  it  was  observed  that  both  exhibited  unstable 
solutions  as  evidenced  by  the  uniform  growth  of  the  error 
in  one  direction  while  that  of  the  Nystrom  and  Euler  had 
decreasing  error  as  the  range  of  iteration  increases.   Be- 
tween the  Euler  and  Nystrom,  the  former  produced  more  accu- 
rate results  though  it  was  a  fraction  of  a  second  longer  in 
computing  time  than  the  Nystrom. 
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TABLE  17-B 

ODE  VI  at  x  =  10.0 

ABSOLUTE  ERROR  FOR  P-C-V,  P-C-VI,  P-C-VII,  and  P-C-VIII 

Double  Precision 


1/h 

EHA 

ZA2 

ZA3 

ZA4 

128 

64 

32 

16 

8 

4 

2 

CPTS 

3.848-10"3 

1.073-10"1 

4.256 

4.556 
-26.070 
-129.471 
-866.001 
5.0 

-48.786 
-99.150 
-191.895 
-347.874 
-557.063 
-613.784 
-52.022 
6.68 

-16.423 
-38.845 
-48.587 
-45.839 
188.851 
295.289 
3621.438 
7.18 

-9.371-10"2 
-9.692-10"1 
-4.337 
-34.281 
-100.571 
-255.149 
-897.303 
4.79 
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TABLE  18 

TRUE  SOLUTION  VALUES  FOR  ODE  VII 
Double  Precision 


v 

J  exact 


1.0  7.615941-10"1 


2.0  9.640275-10"1 


3.0  9.950547-10"1 


4.0  9.993292-10"1 

5.0  9.999092-10"1 


6.0  9.999877,10~1 


7.0  9.999983-10"1 


8.0  9.999997-10"1 

9.0  9.999999-10"1 


10.0  §.^^^'\§~Y 
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TABLE  18 -A 

ODE  VII  at  x  =  10.0 

ABSOLUTE  ERROR  FOR  P-C-I,  P-C-II,  P-C-III,  and  P-C-IV 

Double  Precision 


1/h 


'M 


"N 


'HE 


128 

64 

32 

16 

8 

4 

2 

CPTSl 


5.91-10 


-2 


1.18551-10 


-9 


2.092-10 


11 


6.824836*10 


-1.56-10 


6.78*10 


12 


12 


7.349-10 


-11 


-5.977050-10 


-6 


3.051*10 


-11 


2.5577-10 


8.8764-10 


10 


10 


3.58311-10 


-1.07753-10 


4.42 


-9 


6.000472-10 


-3.322982-10 


3.314797'10 


1.429417 


4.75 


-1.4540-10 


7.3452-10 


10 


10 


3.46305*10 


5.593088*10 


1.750057-10 


2.691970*10 


-5 


-5 


-7 


4.14 


3.963316-10 


-1.511858-10 


9.699544*10 


6.187207-10 


3.852296*10 


6.63 
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TABLE  18-B 

ODE  at  h  =  0.5 

ABSOLUTE  ERROR  FOR  P-C-I,  P-C-II,  P-C-III,  P-C-IV 

Double  Precision 


'M 


'N 


'HE 


1.0 
2.0 
3.0 
4.0 
5.0 
6.0 
7.0 
8.0 
9.0 
10.0 


1.153-10 


8.578*10 


2.967-10 


7.349-10 


1.589-10 


-4 


3.208-10 


6.221*10 


9.152-10 


1.828*10 


1.077-10 


1.007-10 


3.266*10 


-3 


-3 


1. 193*10 


2.764*10 


6.232'iO""2 


-1.250*10 
-1.906-10 
-3.584-10 
1.429 


-1 


-1 


-4.193-10 


1.076*10 


6.486-10 


-3 


-1.726*10 


2.717*10 


-4 


-2.363*10 


1.228*10 


-1.132'IQ 


9.497*10 


-5.593*10 


-7 


-7 


2.027-10 


6.422-10 


6.869-10 


-2.874-10 


-4 


-3 


8.340-10 


2.261-10 


-5.534*10 


-2 


1.000*10 


3.718'IQ 


-1 


3.852-10 
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Results  from  Table  18-C  showed  that  the  last  four  methods 
produced  valid  numerical  solutions  with  Hamming  and  Fourth 

Order  Adams  competing  for  best  accuracy.   As  h  increases 

-  7         -5 
from  h  =  2    to  h  =  2    the  Hamming  method  is  most  precise 

but  from  h  =  2    to  h  =  2    the  Fourth  Order  Adams  method 

showed  better  accuracy.   Computing  time  for  the  fourth  Order 

Adams  method  is  less  than  the  Hamming  method.   The  Euler 

method   compared  with   these    t\\ro   methods    as    h    increases, 

-  7         -4 
started  at  lesser  accuracy  from  h  =  2    to  h  =  2    but 

maintained  its  accurate  results  up  to  h  =  2    and  yielded 

a  much  more  accurate  numerical  solution  at  h  =  2   .   It  also 

had  the  least  computing  time.   Thus  for  classes  of  ODE 

belonging  to  ODE  VII,  the  recommended  order  of  choice  of 

methods  is  as  follows:   Euler,  Nystrom,  Fourth  Order  Adams, 

Hamming,  Third  Order  Adams,  and  lastly  the  Second  Order 

Adams.   The  Milne  and  Hermite  are  not  considered  and  should 

not  be  used  for  this  particular  type  of  ODE  since  they  are 

both  unstable. 

H.   ODE  VIII 

Table  19  presents  the  true  solution  values  for  ODE  VIII 
from  x  =  1.0  to  x  =  10.0.   Data  analysis  from  Tables  19-A 
and  19-B  indicates  that  any  of  the  eight  P-C  sets  can  be 
used  to  solve  ODE  VIII  if  no  specific  criterion  for  accuracy 
is  needed,  since  each  method  yields,  good  and  valid  numerical 
solutions.   However,  for  purposes  of  comparison,  the  Hamming 
method  clearly  stands  out  to  be  the  best  in  terms  of  accuracy 
with  Milne  coming  in  second.   The  Milne  method  provides  the 
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TABLE  18- C 

ODE  at  x  =  10.0 

ABSOLUTE  ERROR  FOR  P-C-V,  P-C-VI,  P-C-VII,  and  P-C-VIII 

Double  Precision 


1/h 


'HA 


'A2 


'A3 


'A4 


128 
64 
32 
16 


4 
2 
CRTS 


0 

0 
5.665862-10 
5.374510-10 
2.632371-10 
1.049942-10 
5.15 


-2.398169-10 


1.929706-10 


-2 


-7 


-6 


1.561650-10 


1.277970-10 


1.068471-10 


9.306166*10 


-3 


-3 


8. 771719*10 


7.05 


-5.968010" 10 
-1.07296-10"9 
-1.87211-10'9 
-3.04922-10"9 
-1.119715*10" 
-3.729009*10" 
-2.033756*10" 
6.20 


10 


-5.0-10 


14 


1.36-10 


-12 


1.715985-10 


-7 


1.727044-10 


-9.570674-10 


4.82 
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TABLE  19 

TRUE  SOLUTION  VALUES  FOR  ODE  VIII 
Double  Precision 


^exact 


1.0  2.319776 

2.0  2.482577 

3.0  1.151562 

4.0  4.691641*10 


1 
-1 


5.0  3.833049*10 

6.0  '    7.562256-10"1 

7.0  .    1.928970 
8.0  2.689507 

9.0  1.510013 


10.0  5.804096-10"1 
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TABLE  19 -A 

ODE  at  x  =  10.0 

ABSOLUTE  ERROR  FOR  P-C-I,  P-C-II,  P-C-III,  and  P-C-VI 

Double  Precision 


1/h 


'M 


'N 


'HE 


128 

64 

32 

16 

8 

4 

2 

CPTS 


1.391257-10 


-9.332-10 


11 


6.532096-10 


-6 


1.850381-10 


7.305628*10 


2.827282-10 


9.532272-10 


-4 


1.051152*10 


6.14 


2.17509-10 


-5.498490-10 


-6.440709-10 


4.145628*10 


9.793708-10 


2.707671*10 


4.63 


7.611170-10 


5.204737-10 


1.836626-10 


-5 


1.595845-10 


5.639547*10 


4.499001*10 


7.254094*10 


2.733740-10 


-4 


3.580536-10 


2.841152-10 


-5 


-4 


7.115245*10 


1.087262*10 


2.270258-10 


5.11 


1.866206-10 


6.48 


-2 


140 


TABLE  19-B 

ODE  at  x  =  10.0 

ABSOLUTE  ERROR  FOR  P-C-V,  P-C-VI,  P-C-VII,  and  P-C-VIII 

Double  Precision 


1/h 

ZHA 

V 

Lk2 

EA3 

ZA4 

128 

64 

32 

16 

8 

4 

2 

CPTS 

- 1  3 
-9.7-10  10 

-4.057-10"11 

-1.838-10"9 

-1.006735-10"7 

-1.414250-10"6 

-7.603359-10"5 

-9.956171-10"4 

5.31 

9.902758-10"7 

3.303524-10"6 

7.843814-10"6 

-1.359301-10"5 

-4.482144-10"4 

-6  .  125271 • 10"2 

-6 . 125271  * 10"2 

7.61 

2.367453-10"7 
1.886624-10"6 
1.498027-10"5 
1.181971'10"4 
9.248289-10"4 
5.783513-10"2 
5. 783513* 10" 2 
8.11 

-2.8983'10"10 
-5.53331-10"9 
-1.155997-10"7 
-1.563035-10"6 
-2.197783'10"5 
5.113170-10"3 
5.113170,10"3 
5.09 
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least  computing  time  with  Hamming  second.   Thus  for  problems 
of  the  type  ODE  VIII,  the  Hamming  method  is  the  most  highly 
recommended,  followed  by  the  Milne,  Fourth  Order  Adams, 
Euler,  Hermite,  Nystrom,  Third  Order  Adams,  and  lastly  the 
Second  Order  Adams. 

I.   ODE  IX 

True  solution  values  are  listed  in  Table  20  for  ODE  IX 
for  the  range  of  integration  x  =  1.0  to  x  =  10.0.   From 
Tables  20-A  and  20-B  it  was  observed  by  analyzing  the 
results  that  all  eight  methods  provided  valid  numerical  solu 
tions.   In  terms  of  better  accuracy  on  the  average  and  the 
least  computing  time  the  Milne  method  is  the  best  candidate 
with  Hamming  coming  in  next,  followed  by  the  Fourth  Order 
Adams,  Hermite,  Third  Order  Adams,  Euler,  Nystrom  and 
lastly  the  Second  Order  Adams.   Thus  if  choice  is  to  be 
made  for  solving  classes  of  ODE  belonging  to  ODE  IX,  the 
Milne  method  is  highly  recommended  with  Hamming  as  the  next 
alternate . 

J.   ODE  X 

This  differential  equation  is  an  example  of  a  controlled 
variable  problem,  where  one  variable  is  expressed  as  a  func- 
tion of  x  and  substituted  into  the  differential  equation  of 
the  other  variable  to  form  a  single  first  order  ODE.   The 
original  two  variable  equations  are: 

y^  =  1.38y1  -  0.81y2  (2-74) 

y£  =  2.16y1  -  1.92y2  (2-75) 
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TABLE  2  0 

TRUE  SOLUTION  VALUES  FOR  ODE  IX 
Double  Precision 


y 


exact 


1.0  2.069535 

2.0  2.249705 

3.0  4.179309 

4.0  9.462527 

5.0  10.633343 

6.0  17.564095 

7.0  42.421352 

8.0  50.806493 

9.0  74.608406 

10.0  186.463649 


14? 


TABLE  20-A 

ODE  IX  at  x  =  10.0 

ABSOLUTE  ERROR  FOR  P-C-I,  P-C-II,  P-C-III,  and  P-C-IV 

Double  Precision 


1/h 

h 

SM 

ZN 

ZHE 

128 

64 

32 

16 

8 

4 

2 

CPTS 

-1.591-10"3 
-6.435*10~3 
-2.578-10"2 
-1.025-10"1 
-4.022-1-"1 
-1.450 
-3.737 
6.99 

1.791-10"8 
1.666-10"7 
2.341-10"6 
3.806'10"5 
1.016-10"3 
4.234-10"2 
9.791-10"1 
5.40 

-1.658-10"3 
-6.522-10"3 
-2.636-10"2 
-1.079-10"1 
-4.569-10"1 
-2.090 
-10.452 
7.09 

1.403-10"5 
1.173-10"4 
9.226-10"4 
7.170-10"3 
5.390-10"2 
3.761-10"1 
2.559 
6.42 
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TABLE  2  0 -B 

ODE  IX  at  x  =  10.0 

ABSOLUTE  ERROR  FOR  P-C-V,  P-C-VI,  P-C-VII,  and  P-C-VIII 

Double  Precision 


1/h 

ZHA 

ZA2 

"A3 

ZA4 

128 

64 

32 

16 

8 

4 

2 

CPTS 

-2.795-10"9 
-6.428-10"8 
-2.136-10"6 
-8.656-10"5 
-9.897-10"4 
-5.696-10"2 
2.108 
5.  73 

-1.648-10"3 
-6.720-10"3 
-2. 788-10"2 
-1.190-10"1 
-5.288* 10"1 
-2.433 
-10.866 
7.93 

6.166-10"5 
4.827-10"4 
3.693-10"3 
2.690-10"2 
1.755-10"1 
8.959-10"1 
3.597 
8.33 

6.922-10"8 
9.598-10"7 
1.610-10"5 
2.948-10"4 
7.000-10"3 
2.261-10"1 
4.380 
5.44 
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with  initial  conditions 

y1(0)  =  -2.9995  (2-76) 

y2(0)  =  4.0010  (2-77) 

and  the  analytic  solutions  given  by 

y1  =  0.0005  e~3x  -  3e~0,3x  (2-78) 

y2  =  0.001  e"3x  +  4e~0,3x  (2-79) 

Since  the  analytical  solutions  of  both  differential  equations 
are  known,  it  is  possible  to  treat  one  variable  as  a  function 
of  x  and  substitute  it  to  the  differential  equation  of  the 
other  variable  to  form  a  single  first  order  ODE  which  now  can 
be  solved  by  the  P-C  sets.   Choosing  y~  as  the  known  function 
given  by  (2-79)  and  substituting  it  in  (2-74),  the  resulting 
equation  is  ODE  X  with  initial  condition  given  by  (2-76)  and 
exact  solution  by  (2-78).   Note  that  y?  remains  a  function 
of  the  single  independent  variable  x  and  as  such  changes  as 
x  changes.   The  other  variable  y,  forms  a  first  order  dif- 
ferential equation  in  the  standard  form  y,'  =  f(x,y).   The 
point  that  is  being  illustrated  here  is  that  this  concept  can 
be  applied  to  problems  like  rate  of  chemical  reaction,  falling 
bodies,  aircraft  or  ballistic  missiles  flight  wherein  time 
variable  is  the  most  important  consideration.   All  other 
variables  can  be  expressed  as  functions  of  the  single  inde- 
pendent variable  time  and  given  constants  or  initial  condi- 
tions.  In  such  cases  the  problem  can  be  reduced  to  a  single 
differential  equation  and  the  methods  considered  herein  can 
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be  applied  (time  (t)  variable  as  x) .   This  concept  does  not 
discount  the  fact  that  these  P-C  methods  can  be  extended  to 
solve  simultaneous  differential  equations  with  some  modifi- 
cations in  the  algorithms.   However,  since  this  extension 
is  not  relevant  to  the  purpose  of  this  particular  paper  in 
analyzing  and  comparing  the  predictor-corrector  methods  as 
applied  to  a  variety  of  ODEs ,  it  is  not  considered  here, 
but  will  be  mentioned  as  a  further  field  of  study. 

Returning  now  to  the  solution  of  ODE  X,  Table  21  shows 
the  true  solution  values  from  x  =  1.0  to  x  =  10.0.   From 
Table  21-A  the  Milne  seemed  to  produce  the  best  accuracy 
but  by  analyzing  further  the  behavior  of  these  errors  in 
step-by-step  integration,  which  is  shown  in  Table  21-B,  it 

i 

is  noted  that  the  errors  of  the  Milne  method  are  increasing 
in  a  uniform  fashion  such  that  as  the  range  of  integration 
is  increased  the  error  grows  while  the  true  solution  values 
as  shown  in  Table  21  decrease.   Thus  eventually  the  error 
will  overcome  the  solution.   The  same  is  true  for  the 
Hermite  method,  while  for  the  Euler  and  the  Nystrom  the 
error  decreases  as  the  range  of  integration  increases  thus 
providing  a  valid  numerical  solution.   The  behavior  of  the 
Milne  and  Hermite  methods  is  expected  since  for  ODE  X  it  is 
easy  to  see  that  C  <  0,  and  as  such  the  Milne  and  Hermite 
methods  are  unstable.   Between  the  Euler  and  the  Nystrom, 
the  Euler  method  is  an  easy  choice  both  in  accuracy  and 
computing  time. 
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TABLE  21 

TRUE  SOLUTION  VALUES  FOR  ODE  X 
Double  Precision 


^exact 


1.0  -2.222429 

2.0  -1.646433 

3.0  -1.219708 


4.0  -9.035826-10 


1 


5.0  -6.693904-10"1 


6.0  -4.958966-10"1 


7.0  -3.673692-10"1 


8.0  -2.721538'10_1 


9.0  -2.016165-10"1 


10.0  -1.493612-10"1 
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TABLE  21 -A 
ODE  at  X  =  10.0 
ABSOLUTE  ERROR  FOR  P-C-I,  P-C-II,  P-C-III,  and  P-C-IV 

Double  Precision 


1/h 


'M 


'N 


'HE 


128 

64 

32 

16 

8 

4 

2 

CPTS 


4.839775*10 


-9.887212*10 


4.893782*10 


1.149738-10 


3.757046-10 


9.800289-10 


1.205*10 


11 


5.3459-10 


-10 


2. 784363*10 


7.03 


-4 


-3.606138-10 


4.286345-10 


1.  576153*10 


2.583136-10 


-8.436503*10 


5.33 


-3.684268-10 


8.584152-10 


3.832304'10"7 


-1.793736-10 


4.903244-10"6 


2.093989-10 


1. 350331-10 


10.0 


-9.818242*10 
-3.913399-10 
-1.558160-10 
-6.236379-10 
-2.594684*10 
-1.288607 
-11.095378 
9.83 
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TABLE  21-B 

ODE  X  with  h  =  0.5 

ABSOLUTE  ERROR  FOR  P-C-I,  P-C-II,  P-C-III,  and  P-C-IV 

Double  Precision 


JM 


JN 


'HE 


1.0 
2.0 
3.0 
4.0 
5.0 
6.0 
7.0 
8.0 
9.0 
10.0 


2.187-10 


2.565-10 


2.144*10 


-3 


1.651*10 


1.239-10 


9.223*10 


6.843'IQ 


5.072-10 


-4 


3.758*10 


2.784*10 


-3.821-10 


-1.905*10 


-3.765-10 


-4 


6.420-10 


-1.076*10 


1.800*10 


-3.012*10 


-3 


-3 


-3 


5.041*10 


8.436*10 


5.056*10 


-1.101*10 


1.016-10 


-4 


-3 


-3 


7.984-10 


6.013-10 


-4 


4.475*10 


3.319-10 


-3 


2.460-10 


1.822*10 


1.350*10 


-5.661-10 

-8.576-10 

-1.542-10 

-2.833-10 

-5.219-10 

-9.619*10 

-1.772 

-3.266 

-6.020 

-11.095 


-2 


-1 


150 


From  Table  21-C,  all  the  last  four  methods  exhibited 
stable  numerical  solutions.   For  accuracy  the  Hamming  method 
offers  the  best  choice,  with  Fourth  Order  Adams  coming 
next,  though  the  computing  time  for  Hamming  is  a  little 
longer  than  the  Fourth  Order  Adams.   Therefore  for  problems 
in  the  class  of  ODE  X,  the  order  of  recommended  preference 
is:   Hamming,  Fourth  Order  Adams,  Third  Order  Adams,  Euler, 
Nystrom,  and  the  Second  Order  Adams  as  the  last  choice.   The 
Milne  and  Hermite  methods  should  not  be  used  for  this  par- 
ticular type  of  problem. 

K .   SUMMARY 

-  7         - 1 
For  series  of  h  values  from  h  =  2    to  h  =  2    and  range 

of  integration  up  to  x  =  10.0,  Tables  22,  22-A,  and  22-B 

list  the  summary  of  results  for  the  eight  predictor-corrector 

sets  considered,  each  run  on  the  test  ODEs  I  to  X.   From 

these  comparative  results,  it  is  clearly  established  that 

the  best  numerical  methods  in  order  of  decreasing  efficiency 

in  general  are: 

1.  Hamming  Modified  P-C  Set 

2.  Fourth  Order  Adams -Moulton  P-C  Set 

3.  Milne  P-C  Set 

4.  Third  Order  Adams -Moulton  P-C  Set 

5.  Euler  P-C  Set 

6.  Nystrom  P-C  Set 

7.  Second  Order  Adams -Moulton  P-C  Set 

8.  Hermite  P-C  Set. 
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TABLE  21-C 

ODE  X  at  x  =  10.0 

ABSOLUTE  ERROR  FOR  P-C-V,  P-C-VI,  P-C-VII,  and  P-C-VIII 

Double  Precision 


1/h 


'HA 


'A2 


'A3 


'A4 


128 
64 
32 
16 


4 
2 
CPTS 


9.0*10 


14 


8.773099-10 


7.065215-10 


4.17081*10 


1.149665-10 


-7 


8.476336*10 


-6.538981-10 


-5.207967*10 


-4.272440*10 


-4 


-4.137219-10 


6.02 


-5 


3.630456-10 


3.165651*10 


9.45 


9.07653-10 


7.319508-10 


5.949594-10 


4.912304*10 


4.180375-10 


-3. 757718*10 


3.628554*10 


9.98 


-4.0-10 


14 


6.9-10 


-13 


1.380-10 


-11 


3.0945*10 


-10 


2.697520-10 


1.553614-10 


1.293084-10 


5.45 


-6 
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TABLE  2  2 

SUMMARY  OF  RESULTS  FOR  ODE  I  TO  ODE  IV 

FOR  h=2~7  TO  h=2_1  AT  x=10.0 


ODE 
Number 

Stable  P-C  Sets  in 
Order  of  Preference 

Average 
Accuracy 

Computing 
Time 
in  sec. 

Unstable 
P-C  Sets 

I 

1. 

Hamming 

IO"4 

5.76 

1.  Milne 

2. 

Fourth  Order  Adams 

IO"3 

5.99 

2.  Hermite 

3. 

Third  Order  Adams 

IO"3 

8.55 

4. 

Euler 

IO"2 

7.65 

5. 

My strom 

IO"2 

7.46 

6. 

Second  Order  Adams 

IO-2 

7.98 

II 

1. 

Hamming 

10"4 

5.26 

1.  Hermite 

2. 

Fourth  Order  Adams 

10"3 

4.96 

3. 

Milne 

-  3 

10  ° 

4.69 

4. 

Euler 

io-1 

6.38 

5. 

Second  Order  Adams 

10"1 

7.79 

6. 

Third  Order  Adams 

10"1 

8.36 

7. 

My  strom 

10"1 

5.34 

III 

1. 

Hamming 

IO"8 

6.23 

None 

2. 

Milne 

IO"7 

5.59 

3. 

Fourth  Order  Adams 

10"7 

8.23 

4. 

Third  Order  Adams 

IO-5 

10.0 

5. 

Hermite 

10"4 

7.92 

6. 

Euler 

10"4 

7.81 

7. 

Second  Order  Adams 

10"4 

10.0 

8. 

Nystrom 

10"4 

7.12 

IV 

1. 

Hamming 

IO"10 

4.45 

None 

2. 

Milne 

10-9 

3.83 

3. 
4. 

Fourth  Order  Adams 
Third  Order  Adams 

io-8 

IO"6 

4.02 
7.02 

5. 

My  strom 

io"5 

3.67 

6. 

Second  Order  Adams 

IO"5 

6.46 

7. 

Euler 

IO"5 

4.39 

8. 

Hermite 

IO"4 

6.74 
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TABLE  22-A 

SUMMARY  OF  RESULTS  FOR  ODE  V  TO  ODE  VIII 

FOR  h=2~7  TO  h=2_1  AT  x=10.0 


ODE 

Stable  P-C  Sets  in 

Computing 
Time 

1 
1 

i 

1 

Unstable 

Number 

0 

rder  of  Preference 

Accuracy 

j 

P-C  Sets 

in  sec. 

1 

1 
1 

V 

1. 

Hamming 

10-13 

6.31 

None 

2. 

Milne 

10"12 

5.58 

3. 

Fourth  Order  Adams 

10-12 

7.78 

4. 

Third  Order  Adams 

io-n 

10.0 

5. 

Hermite 

10-1° 

5.35 

6. 

^ystrom 

10-7 

5.16 

7. 

Euler 

10-7 

5.63 

8. 

Second  Order  Adams 

10"7 

10.0 

VI 

1. 

Hamming 

lO"2 

5.0 

1. 

Euler 

2. 

Fourth  Order  Adams 

10"! 

4.79 

2. 

Nys trom 

3. 

Milne 

10"! 

4.40 

3. 
4. 

5. 

Hermite 

Third  Order 
Adams 

Second  Order 
Adams 

VII 

1. 

Euler 

lO"" 

4.42 

1. 

Milne 

2. 

Nys trom 

-9 

10  y 

4.14 

2. 

Hermite 

3. 

Fourth  Order  Adams 

10"8 

4.82 

4. 

Hamming 

10"8 

5.15 

5. 

Third  Order  Adams 

10-7 

6.20 

6. 

Second  Order  Adams 

10-4 

7.05 

VIII 

1. 

Hamming 

10'8 

5.31 

None 

2. 

Milne 

10-7 

4.63 

3. 

Fourth  Order  Adams 

10"6 

5.09 

4. 

Euler 

10"5 

6.14 

5. 

Hermite 

10"5 

6.48 

6. 

Nys trom 

10-4 

5.11 

7. 

Third  Order  Adams 

10-4 

8.11 

8. 

Second  Order  Adams 

10-4 

7.61 
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TABLE  22 -B 
SUMMARY  OF  RESULTS  FOR  ODE  IX  AND  ODE  X 
FOR  h=2~7  TO  h=2_1  AT  x=10.0 




ODE 

Number 

Stable  P-C  Sets  in 
Order  of  Preference 

Average 

Accuracy 

Computing 
Time 
in  sees . 

Unstable 
P-C  Sets 

IX 

1.  Milne 

lO"4 

5.40 

2.  Hamming 

10"4 

5.73 

3.  Fourth  Order  Adams 

10"4 

5.44 

4.  Hermite 

10"3 

6.42 

5.  Third  Order  Adams 

10 

8.33 

6.  Euler 

i(T2 

6.99 

7.  'Ny strom 

10"2 

7.09 

8.  Second  Order  Adams 

10"2 

7.93 

X 

1.  Hamming 

lO"10 

6.02 

1.  Milne 

2.  Fourth  Order  Adams 

-  9 
10 

5.45 

2.  Hermite 

3.  Third  Order  Adams 

i(T6 

9.98 

4.  Euler 

10"6 

7.03 

5.  Nystrom 

10_6 

10.0 

6.  Second  Order  Adams 

10"5 

9.45 
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The  most  important  criteria  applied  is  accuracy  followed  by 
computing  time.   It  is  obvious  that  the  Hamming  method  did 
not  produce  the  least  computing  time  as  compared  to  the 
other  methods  especially  with  that  of  Milne  and  the  Fourth 
Order  Adams  methods.   This  difference  of  less  than  a  second 
in  computing  time  is  attributed  to  the  fact  that  the  Hamming 
method  used  additional  computation  for  the  truncation  error 
to  modify  the  predicted  and  corrected  values  for  every 
iteration.   Analyzing  the  computing  time  for  each  method 
within  the  range  of  their  stability,  the  list  below  shows 
the  order  of  increasing  computing  time: 

Average 
P-C  Set Computing  Time 

1.  Milne  4 . 82  sees . 

2.  Hamming  5.60  sees. 

3.  Nystrom  5.65  sees. 

4.  Fourth  Order  Adams  6.11  sees. 

5.  Euler  6.19  sees. 

6.  Hermite  6.58  sees. 

7.  Second  Order  Adams  8.40  sees. 

8.  Third  Order  Adams  8.69  sees. 

This  confirmed  the  idea  that  convergence  of  a  corrector  does 

not  necessarily  assure  the  convergence  towards  the  true 

solution  values,  y(x  ) ,  but  only  to  some  definite  value, 

y  ^, .   This  is  best  illustrated  in  the  case  of  Hermite  and 
'  n+1 

Third  Order  Adams  methods  with  6.58  sees  and  8.69  sees  average 
computing  time,  respectively.   By  recalling  from  the  order 
of  efficiency  of  the  P-C  sets,  the  Third  Order  Adams  ranks 
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fourth  in  accuracy  while  Hermite  ranks  eighth,  which  points 
to  the  fact  that  though  the  Hermite  method  converges  much 
faster  than  the  third  order  Adams,  its  accuracy  is  much 
poorer  than  the  Third  Order  Adams.   Clearly  then  it  can  be 
seen  that  the  Hermite  method  is  converging  only  to  some 
definite  value,  y    ,  but  not  to  the  true  solution,  y(*  )> 
otherwise  it  should  have  much,  better  accuracy  than  the  Third 
Order  Adams  method,  which-  took.  much,  longer  to  converge  to 

a  definite  value  y  L, ,  in  each  iteration.   This  analysis  is 

' n+1 7  ' 

meaningful  due  to  the  fact  that  a  standard  mode  of  applica- 
tion  is  used,  that  of  P(EC)  .   It  should  be  noted  however 
that  rapid  convergence  is  essential  to  accuracy  as  shown 
by  the  Hamming,  Milne,  and  Fourth  Order  Adams  methods  which 
all  rank  well  above  the  other  methods  in  both  accuracy  and 
computing  time  in  general. 
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X.   CONCLUSIONS 

The  following  conclusions  can  be  drawn  and  recommenda- 
tions made  from  the  analysis  and  numerical  results  obtained 
in  this  paper. 

1.  Numerical  instability  results  from  extraneous  solutions 
of  the  difference  equations  which  bear  no  connection   to 
the  exact  solution.   The  conditions  of  asymptotic  (strong) 
and  absolute  stability  can  be  seen  clearly  by   reference 
to  the  extraneous  solutions.   If  in  the  limit  h  -*■  0 ,  the 
extraneous  solutions  vanish  as  n  ->  °°  ,  then  the  method  is 
strongly  stable  and  convergent.   If,  for  values  of  h  less 

than  some  h  ,  the  extraneous  solutions  vanish  as  n  -*■  «>  the 

o'  ' 

method  is  absolutely  stable.  Relative  stability  is  a  sig- 
nificant concept  for  ODE  where  f  >  0.  The  condition  pro- 
vides that  extraneous  solutions  will  not  grow  more  rapidly 
or  decay  more  slowly  than  the  true  solution.  Thus,  before 
a  method  is  used,  the  characteristic  roots  of  its  related 
difference  equation  must  be  analyzed. 

2.  The  finite  real  negative  stability  bounds  for  the  P-C 
methods  considered,  as  determined  by  numerical  experiments, 
are  in  reasonable  agreement  with  theoretical  predictions. 
As  such  experimental  bounds  could  be  used  as  a  guide  to  the 
proper  selection  of  a  method  to  solve  a  particular  problem. 

3.  The  stability  of  a  P-C  set  depends  on  both  the  predictor 
and  corrector  equations.   Thus,  when  two  P-C  sets  have  the 
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same  corrector  but  different  predictors,  choose  the  one  with 
a  better  predictor.   Likewise,  when  two  P-C  sets  use  the 
same  predictor  with  different  correctors,  choose  the  one 
with  a  better  corrector. 

4.  P-C  sets  with  higher  order  truncation  errors  are  not 
necessarily  better.  In  choosing  a  P-C  set,  the  order  of 
truncation  error  should  not  be  the  sole  criterion. 

5.  The  convergence  of  the  corrector  formula  will  ensure 
that  the  sequence  of  approximations  will  converge  to  some 
definite  value  but  not  necessarily  to  the  true  solution 
values.   As  such,  the  fast  computation  time  of  a  method 
does  not  necessarily  imply  greater  accuracy  for  the  result- 
ing numerical  solution. 

6.  If  the  integration  will  involve  a  large  number  of  steps, 
a  stable  method  should  be  used. 

7.  If  the  function  evaluation  is  lengthy,  the  range  of  inte- 
gration large,  and  better  accuracy  and  lesser  cost  of  com- 
puting time  are  prime  considerations,  then  the  following  P-C 
methods  are  recommended,  based  on  the  overall  efficiency 
they  have  demonstrated  in  numerical  experiments  on  a  wide 
variety  of  different  test  ODEs : 

a)  Hamming  P-C  Sets 

b)  Fourth  Order  Adams  P-C  Set 

c)  Milne  P-C  Set 

d)  Third  Order  Adams  P-C  Set 

e)  Euler  P-C  Set 

f)  Nystrom  P-C  Set 
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g)   Second  Order  Adams  P-C  Set 
h)   Hermite  P-C  Set 

The  listing  indicates  order  of  preference.   However  caution 
must  be  exercised  not  to  use  Milne  and  Hermite  methods  if 
fy  <  0. 

8.  For  maximum  accuracy  to  be  achieved,  select  a  step  size 
h,  for  which  the  truncation  error  and  the  roundoff  error  have 
equal  orders  of  magnitude. 

9.  If  high  accuracy  results  with  less  roundoff  error  are 
desired,  the  procedure  should  use  double  precision,  which 
involves  only  slightly  more  computing  time  than  single  pre- 
cision in  an  IBM  360/67  computer. 
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XI.   EXTENSIONS 

The  following  topics,  which  are  direct  extensions  of 
this  paper,  are  interesting  areas  of  further  research  and 
study  in  the  wide  field  of  numerical  analysis: 

1.  Stability  versus  Accuracy 

A  necessary  condition  for  numerical  integration  is 
stability*  that  is,  a  stable  value  of  h  must  be  employed 
when  a  fhite  stability  boundary  exists.   Nevertheless, 
there  is  no  guarantee  that  all  values  of  h  from  zero  to  the 
limiting  \aLue  will  yield  accurate   results.      Thus  the 
question  that  must  be  faced  is:   Of  what  value  is  a  method 
that  is  stable  in  a  region  where  it  is  inaccurate?   The 
answer  to  this  question  will  show  that  the  relationship 
between  stability  and  accuracy  is  fundamental  to  the  choice 
of  h  for  a  particular  method. 

2.  Increasing  the  Stability  Bounds 

It  has  been  shown  that  in  some  methods  the  real  stability 
bound  is  somewhat  constrained.   Two  possible  ways  where  the 
stability  bounds  can  be  increased  are: 

a)  By  the  use  of  a  weighting  factor  which  involves  experi- 
menting with  different  values  of  the  free  parameters  that 
are  used  in  the  method  by  undetermined  coefficients,  similar 
to  Hamming's  corrector  formula  derivation. 

b)  By  averaging,  which  means  computing  the  new  value  of 
yn+,  after  a  finite  number  of  steps  (say  fifty)  of  a  particular 
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method,  as  the  average  of  the  old  y  _,,  and  a  value  called 
'  to  J  n+1 

y  +,  generated  by  another  method.   This  should  be  done 
periodically. 

3.  Predictor-Corrector  Methods  Interaction 

The  interaction  of  the  different  methods  in  solving  a 
particular  problem  involves  use  of  the  P-C  sets  as  sub- 
routines, wherein  one  method  will  be  used  to  provide  the 
solution  for  smaller  values  of  h,  then  another  method  will 
be  used  to  solve  for  larger  values  of  h.   This  should  be 
quite  interesting  since  it  has  been  shown  that  some  methods 
exhibit  better  accuracy  for  small  values  of  h,  then  become 
inaccurate  for  large  values  of  h,  while  other  methods  main- 
tain their  accuracy  up  to  large  values  of  h  but  might  be 
prohibitive  to  use  for  small  values  of  h  due  to  their  longer 
computing  time. 

4.  Solution  of  Systems  of  Differential  Equations 
Revision  of  the  algorithms  to  enable  the  solution  of 

systems  of  differential  equations  and  higher  order  differ- 
ential equations.   This  extension  is  straightforward.   It 
involves  only  extra  computational  steps  using  the  same 
formulas.   Complexity  arises  from  the  need  to  use  vectors 
and  arrays  for  temporary  storage. 

5.  Practical  Applications 

Applications  to  real-life  problems  such  as  heat  flow 
problems,  simple  electrical  circuits,  force  problems,  rate 
of  bacterial  growth,  rate  of  decomposition  of  radioactive 
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material,  crystallization  rate  of  a  chemical  compound,  rate 
of  population  growth,  and  so  forth. 

6.   Adjustment  of  the  Step  Size  During  the  P-C  Solution 

This  involves  monitoring  the  modified  truncation  error 
developed  by  Hamming.   If  the  value  is  too  large,  then 
the  step  size  is  too  large,  and  the  calculation  should  be 
repeated  with  a  smaller  value  of  h,  say  h/2.   Note  that  it 
is  not  necessary  to  go  back  to  the  beginning  of  the  calcula 
tion  but  only  to  the  point  at  which  the  truncation  error 
becomes  too  large.   If  the  truncation  error  is  too  small, 
computation  time  is  being  wasted  and  h  should  be  increased, 
say  to  2h. 
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FLOW  CHART  FOR  ALGORITHM  1 


INPUT  TERMNNAL  BOUNDARY, 
NUMBER  OF  ITERATIONS, 
CONVERGENCE  TERM,  NUMBER 
OF  STARTING  POINTS 


& 


READ  STEP  SIZE  AND 
CONTROL  PRINTING 


SPECIFY  INITIAL  CONDITION 


Yes 


-J    STOP 


PRINT  INPUT  DATA  SET 


COMPUTE  FUNCTIONAL 
EVALUATIONS  OF 
STARTING  POINTS 


CALL  FOURTH  ORDER 

RUNGE-KUTTA  METHOD 

TO  GENERATE 

NEEDED  ADDITIONAL 

STARTING  POINTS 
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COMPUTE  PREDICTED 
SOLUTION  VALUE  USING 
PREDICTOR  FORMULA 


B 


COMPUTE  FUNCTIONAL 
EVALUATION  OF  PREDICTED 
SOLUTION  VALUE 


COMPUTE  CORRECTED 
SOLUTION  VALUE  USING 
CORRECTOR  FORMULA 


COMPUTE  MAGNITUDE 

OF  DIFFERENCE  BETWEEN 

PREDICTED  VALUE  AND 

CORRECTED  VALUE 


SET  PREDICTED 
SOLUTION  VALUE 
EQUAL  TO  CORRECTED 
SOLUTION  VALUE 


B 
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COMPUTE  TRUE 
SOLUTION  VALUE 


COMPUTE  NUMERICAL 
ERROR  IN  CALCULATION 


JL 


PRINT 
RESULTS 


UPDATE  REQUIRED 
PARAMETER  VALUES 


I 


ADVANCE  INTEGRATION 
POINT  BY  A  STEP 
SIZE  INCREMENT 


Yes 


No 


.0 
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FLOW  CHART  FOR  ALGORITHM  2 


INPUT  TERMINAL  BOUNDARY, 
NUMBER  OF  ITERATIONS, 
CONVERGENCE  TERM,  NUMBER 
OF  STARTING  POINTS 


D 


READ  STEP  SIZE 
AND  CONTROL  PRINTING 


SPECIFY  INITIAL  CONDITION 


Yes 


PRINT  INPUT  DATA  SET 


COMPUTE  FUNCTIONAL 
EVALUATIONS  OF 
STARTING  POINTS 


CALL  FOURTH  ORDER 

RUNGE-KUTTA  METHOD 

TO  GENERATE 

NEEDED  ADDITIONAL 

STARTING  POINTS 
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COMPUTE  PREDICTED 
SOLUTION  VALUE  USING 
PREDICTOR  FORMULA 


COMPUTE  MODIFIED 
PREDICTED  SOLUTION 
VALUE 


COMPUTE  FUNCTIONAL 
EVALUATION  OF  MODIFIED 
PREDICTED  VALUE 


COMPUTE  CORRECTED  SOLUTION 

VALUE  USING  CORRECTOR 
FORMULA  THEN  COMPUTE  MODI- 
FIED CORRECTED  SOLUTION  VALUE 


COMPUTE  MAGNITUDE  OF 
DIFFERENCE  BETWEEN  MODI- 
FIED PREDICTED  VALUE  AND 
MODIFIED  CORRECTED  VALUE 


168 


G> 


COMPUTE  TRUE 
SOLUTION  VALUE 


COMPUTE  NUMERIAAL 
ERROR  IN  CALCULATION 


PRINT 
RESULTS 


UPDATE  REQUIRED 
PARAMETER  VALUES 


ADVANCE  INTEGRATION 
POINT  BY  A  STEP  SIZE 
INCREMENT 


No 


SET  MODIFIED 
PREDICTED  VALUE 
EQUAL  TO  MODIFIED 
CORRECTED  VALUE 


Yesy 


D 


169 


PROGRAM  ONE 

EULER  PREDICTOR-CORRECTOR  METHOD.  ODE:  Y'=-Y 
WITH  Y(0)=1.  TO  SOLVE  ANOTHER  ODE  SIMPLY  CHANGE 
INITIAL  CONDITIONS  AND  FUNCTION  SUBROUTINES. 


C  *  FORTRAN  PROGRAM  FOR  EULER  PREDICTOR-CORRECTOR  * 

C  *  METHOD  IN  THE  NUMERICAL  SOLUTION  OF  THE  ORDINARY  * 

C  *  DIFFERENTIAL  EQUATION  IN  THE  FORM  OF  DY/DX=F(X,Y)   * 

C  *  WITH  THE  INITIAL  CONDITION  Y(XO)=YO  * 

C  *  THE  PARAMETERS  TO  THE  PROGRAM  HAVE  THE  FOLLOWING  * 

C  *  MEANINGS...  * 

C  *    XMAX THE  TERMINAL  BOUNDARY  CONDITION  * 

C  *    EPS THE  CONVERGENCE  TEST  CONSTANT  * 

C  *    MAX  THE  MAXIMUM  NUMBER  OF  ITERATIONS  * 

C  *    YP_Z THE  PREDICTED  VALUE  OF  THE  NEXT  * 

C  *              SOLUTION  POINT  * 

C  *    YC THE  CORRECTED  VALUE  OF  THE  NEXT  * 

C  *              SOLUTION  POINT  * 

C  *    FP THE  PREDICTED  FUNCTION  EVALUATION  OF  YP* 

C  *    FC THE  CORRECTED  FUNCTION  EVALUATION  OF  YC* 

C  *    H ITHE  STEP  SIZE  TO  BE  USED  TO  ADVANCE  * 

C  *              THE  POINT  OF  INTEGRATION  * 

C  *    ICON THE  DESIRED  POINT  OF  PRINTING  THE  * 

C  *              COMPUTED  SOLUTION  FOR  EVERY  FIXED  * 

C  *              NUMBER  OF  INTEGRATION  POINTS  * 

C  *    YEXACT THE  TRUE  SOLUTION  VALUES  * 

C  *    ERROR THE  ABSOLUTE  DEVIATION  OF  THE  NUMERICAL* 

C  *          ""    SOLUTION  FROM  THE  TRUE  SOLUTION  * 

C  *    FCT FUNCTION  SUBROUTINE  TO  EVALUATE  ODE  * 

C  *              FUNCTION  * 

C  *    EXACT FUNCTION  SUBROUTINE  TO  COMPUTE  TRUE  * 

C  *              SOLUTION  VALUES  * 

C 

c 

C      INITIALIZE  INPUT  CONSTANTS 
C 

MAX=2 

XMAX=10.0 

EPS-0. 000001 

c 

C      READ  STEP  SIZE  AND  CONTROL  PRINTING 
C 

5     READ(5,100)  H,ICON 
C 

C      SPECIFY  INITIAL  CONDITION 
C 

XO=0.0 

Y0=1.0 

c 

C      TEST  FOR  END  OF  INPUT  DATA 
C 

IF(ICON-O)  45,45,7 
C 

C      PRINT  INPUT  DATA  SET 
C 
7     WRITE(6,  1000)  H,ICON 

WRITE(6,1001 )  X0,Y0 
C 

C      COMPUTE  NEXT  SOLUTION  POINT 
C 

NC=( XMAX+H/2.) /H 

FP=FCT(X0,Y0) 

DO  40  N=1,NC 
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M=l 

XO=X0+H 

YP=Y0+H*FP 
11    FC=FCT(XO,YP) 

YC=Y0+H/2.*(FP+FC) 

OELY=ABS(YC-YP) 
C 

C      TEST  FOR  CONVERGENCE 
C 

IF(DELY-EPS)  30,30,15 
C 

C  TEST    IF    MAXIMUM    NUMBER    OF    ITERATIONS    IS    SATISFIED 

C 
15  IF(M-MAX)     20,20,30 

20         YP=YC 

M=M+1 

GO    TO    11 
C 

C      COMPUTE  TRUE  SOLUTION  VALUE 
C 

30    YEXACT=EXACT(XO) 
C 

C      COMPUTE  ERROR  IN  CALCULATION 
C 

ERROR=YEXACT-YC 
C 

C      TEST  IF  DESIRED  POINT  CF  PRINTING  IS  REACHED 
C 

IF{ (N/I CON* ICON J.EQ.N)  WRI T E (6 , 100  2 )X0 , YP, YC , Y EXACT , 
1ERR0R 
C 

C      UPDATE  REQUIRED  PARAMETER  VALUES 
C 

YO=YC 

FP=FC 
40         CONTINUE 
C 

C  LOOP    BACK    TO    READ    NEXT.  SET    OF    INPUT    DATA 

C 

GO  TO  5 
45    STOP 
100   FORMAT(F10.7,I10) 

1000  FORMAT C////41X,  '*INPUT  DATA  SET  USED** //43X, «H     =«  , 
1F10.7/43X,' ICON  =',I10) 

1001  FORMAT (//29X, «**RESULTS  OF  EULER  PREDICTCP-COPRECTOR 
1METH0D**«//10X,  '    X   ' ,5X,»    YP    ' ,5X,»    YC    •, 

15X,«    YEXACT    '^X,'    ERROR    •  //12X  ,  F10  .7  ,44X  , 
1F10.7) 

1002  FORMAT ( » 0' ,9X , F15 .7 , 2X ,F15 .7, 2X, Fl 5. 7, 2X, F 15.7, 2X, 
1F15.7) 

END 

FUNCTION  FCT(XN,YN) 

FCT=-YN 

RETURN 

END 

FUNCTION  EXACT(XX) 

EXACT=EXP(-XX) 

RETURN 

END 
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*INPUT  DATA  SET  USED* 
H     =  0.5000000 
ICON  =  2 

**RESULTS  OF  EULER  PREDICTOR-CORRECTOR  METHOD** 


0.0000000 
1.0000000 
2.0000000 
3.0000000 
4.0000000 
5.0000000 
6.0000000 
7.0000000 
8.0000000 
9.0000000 
10.0000000 


YP 

0.3588867 
0.1309212 
0.0477613 
0.0174238 
0.0063564 
0.0023189 
0.0008459 
0.0003086 
0.0001126 
0.0000411 


YC 

0.3634033 
0.1325720 
0.0483635 
0.0176435 
0.0064365 
0.0023481 
0.0008566 
0.0003125 
0.0001140 
0.0000416 


YEXACT 
1.0000000 
0.3678794 
0.1353353 
0.0497871 
0.0183156 
0.0067379 
0.0024788 
0.0009119 
0.0003355 
0.0001234 
0.0000454 


ERROR 

0.0044761 
0.0027633 
0.0014236 
0.0006722 
0.0003014 
0.0001307 
0.0000553 
0.0000230 
0.0000094 
0.0000038 


*INPUT  DATA  SET  USED* 
H     =  1.0000000 
ICON  =  1 

**RESULTS  OF  EULER  PREDICTOR-CORRECTOR  METHOD** 


0.0000000 
1.0000000 
2.0000000 
3.0000000 
4.0000000 
5.0000000 
6.0000000 
7.0000000 
8.0000000 
9.0000000 
10.0000000 


YP 

0.2500000 
0.1562500 
0.0507813 
0.0258789 
0.0095825 
0.0044327 
0.0017519 
0.0007731 
0.0003156 
0.0001361 


YC 

0.3750000 
0.1718750 
0.0683594 
0.0300293 
0.0122986 
0.0052910 
0.0021987 
0.0009362 
0.0003919 
0.0001660 


YEXACT 
1.0000000 
0.3678794 
0.1353353 
0.0497871 
0.0183156 
0.0067379 
0.0024788 
0.0009119 
0.0003355 
0.0001234 
0.0000454 


ERROR 

-0.0071206 
-0.0365397 
-0.0185723 
-0.0117137 
-0.0055606 
-0.0028122 
-0.0012868 
-0.0006007 
-0.0002685 
-0.0001206 
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PROGRAM    TWO 

MILNE    PREDICTOR-CORRECTOR    METHOD.    ODE:    Y»=-Y 
WITH    Y(0)=1.    TO    SOLVE    ANOTHER    ODE    SIMPLY    CHANGE 
INITIAL    CONDITIONS    AND    FUNCTION    SUBROUTINES. 


CV.    -^   J.    -'  -    -J*-    -'-    •Jf   O-   ^1-   V,   ^',   ,',   -V    -L-   -^  ~K.    ^V   *-',    *l,   O-    -J-   ,1-   V-    *•-   S.U  »1  -    »<-   »W    O-   OU    *V    »!-.    U#JU   »',    J,   ^,   J,    .1.    J-    4.    J,   J,   J-   ^.     J.  a   J.   *   J.   Jf   J,   J,  J. 
T  'Is  '.'  '.-  'i*  'i-   'i*  i'  '^  'f  'i*  'i-  t  'i*  T  'i*  ".'  v  i-  -(*  ^  '■'  'c  •>*  -c  V  ',»  •«  Jp  t  T  ^  tt»  i*  *n  *»*  *<*  T  'S-  "V  *t*  T-  1"  <*  t*  *«■•  T*  *?  *r  nr  ***  «T*  *n 

C  *  FORTRAN    PROGRAM    FOR    MILNE    PREDICTOR-CORRECTOR  * 

C  *  METHOD    IN    THE    NUMERICAL    SOLUTION    OF    THE    ORDINARY  * 

C  *  DIFFERENTIAL    EQUATION    IN    THE    FORM    OF    DY/DX=F(X,Y)  * 

C  *  WITH    THE    INITIAL    CONDITION    Y(XO)=YO  * 

C  *  THE    PARAMETERS    TO    THE    PROGRAM    HAVE    THE    FOLLOWING  * 

C  *  MEANINGS...  * 

C  *         XMAX THE    TERMINAL    BOUNDARY    CONDITION  * 

C  *         EPS THE    CONVERGENCE    TEST    CONSTANT  * 

C  *         MAX       ITHE    MAXIMUM    NUMBER    OF    ITERATIONS  * 

C  *         YP THE    PREDICTED    VALUE    OF    THE    NEXT  * 

C  *  SOLUTION    POINT  * 

C  *         YC THE    CORRECTED    VALUE    OF    THE    NEXT  * 

C  *  SOLUTION    POINT  * 

C  *         FP THE    FUNCTION    EVALUATION  * 

C  *         H THE    STEP    SIZE    TO    BE    USED    TO    ADVANCE  * 

C  *  THE    POINT    OF     INTEGRATION  * 

C  *          ICON            THE    DESIRED    POINT    OF    PRINTING    THE  * 

C  *  "*              COMPUTED    SOLUTION    FOR    EVERY    FIXED  * 

C  *  NUMBER    OF    INTEGRATION    POINTS  * 

C  *         YEXACT THE    TRUE    SOLUTION    VALUES  * 

C  *         ERROR THE    ABSOLUTE    DEVIATION    OF    THE    NUMERICAL* 

C  *  SOLUTION    FROM    THE    TRUE    SOLUTION  * 

C  *         FCT FUNCTION    SUBROUTINE    TO    EVALUATE    ODE  * 

C  *  FUNCTION  * 

C  *         EXACT FUNCTION    SUBROUTINE    TO    COMPUTE    TRUE  * 

C  *  SOLUTION    VALUES  * 

C  *         RKUTTA THE    SUBROUTINE    USED    TO    GENERATE  * 

C  *  NEEDED    STARTING    VALUES  * 

C-y-  -v  v>*  **•  V'  -*-*■'-  V'  -**  v-  V'  "■*-  »**  y«»  **»  «.v  y*  •*■  y '  ***  >v  v-  *■■  -  -  -  -»'-  v  -  *■'-  *  ■*  y*  »**  ***  **-  -1*  >  -  *-  -  ■  •*-***  -^  •***■  -*'  v-  -v  -j-  ■**-  -^  v-  v*  ***  Jr  v-*  **r  •**  -'- 
■nr  or  Tr  *v*  n*  *V*  n*  A*1  *»**  -V  *>*  *c*  *¥*  *i*  t-  "ir  -v  •nr  n*  'r*  ^o  n^  -**  -*■*  t*  *r*  n*  -r*  J|»  »*»  -r*  *v*  ^n  'i^  *t*  ******  *i*  *i"*  ***  ^*  ^T"  "r  *i*  *»^  "P  n-  *r  T*  Or  "*P  or  or  or 

C 
C 

C  INITIALIZE    INPUT    CONSTANTS 

C 

XMAX-10.0 

MAX=2 

EPS=0.00C001 

K=3 
C 

C  READ    STEP    SIZE    AND    CONTROL    PRINTING 

C 

5  READ(5,100)    H,ICON 

C 

C  SPECIFY    INITIAL    CONDITION 

C 

XO=0.0 

YO=1.0 

YR  =  Y0 
C 

C  TEST    FOR    END    OF     INPUT    DATA 

C 

IF(ICON-O)    45,45,7 
C 
C  PRINT    INPUT    DATA    SET 


C 

7  WRITE(6,1000)     H,ICON 

WRITE(6,1001)     X0,YO 
C 
C  COMPUTE    NEXT    SOLUTION    POINT 
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NC=(XMAX+H/2. )/H 

X1=X0+H 

X2=X0+2.*H 

X3=X0+3.*H 

X=X0+4.*H 

CALL    RKUTTA(X0,YR,Y1,Y2,Y3,K,H) 

F1=FCT(X1,Y1) 

F2=FCT(X2,Y2) 

F3=FCT(X3,Y3) 

DO    40    N=4,NC 

M=l 

YP=Y0+4.*H/3.*{2-*Fl-F2+2.*F3) 

FP=FCT(X,YP) 
11         YC=Y2+H/3.*(F2+4.*F3+FP) 

DELY=ABS(YC-YP) 
C 

C  TEST    FOR    CONVERGENCE 

C 

IF(DELY-EPS)    30,30,15 
C 

C  TEST    IF    MAXIMUM    NUM8ER    OF    ITERATIONS    IS    SATISFIED 

C 
15         IF(M-MAX)     20,20,30 
20         FP=FCT(X,YC) 

M=M+1 

GO    TO    11 
C 

C      COMPUTE  TRUE  SOLUTION  VALUE 
C 

30    YEXACT=EXACT(X) 
C 

C      COMPUTE  ERROR  IN  CALCULATION 
C 

ERROR=YEXACT-YC 
C 

C      TEST  IF  DESIRED  POINT  OF  PRINTING  IS  REACHED 
C 

IF( (N/ ICON* ICON J.EQ.N).  WRITE (6, 1002)  X, YP , YC , YEXACT, 
1ERR0R 
C 

C      UPDATE  REQUIRED  PARAMETER  VALUES 
C 

F1  =  F2 

F2  =  F3 

F3  =  FP 

Y0  =  Y1 

Y1=Y2 

Y2=Y3 

Y3=YC 

Y—  Y  j.  U 

40    CONTINUE 
C 

C      LOOP  BACK  TO  READ  NEXT  SET  OF  INPUT  DATA 
C 

GO  TO  5 

45    STOP 

100   FORMAT(F10.7,I10) 

1000  F0RMAT(////41X,'*INPUT  DATA  SET  US  ED*' //43X , • H     =• 
1F10.7/43X, f ICON  =■ ,110) 

1001  FQRMAT(//29X,»**RESULTS  OF  MILNE  PREDICTOR-CORRECTOR 
lMETH0D**»//10Xt «    X    '^X,1    YP    ',5X,'    YC    ' 
15X,1    YEXACT    ',5X,«    ERROR    • //12X, F10. 7 ,44X , 
1F10.7) 

1002  FORMAT ( '0',9X,F15.7,2X,F15.7,2X,F15.7,2X,F15.7,2X, 
1F15.7) 

END 
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SUBROUTINE  RKUTT A{ XX , YY, Yl , Y2, Y3 ,KK, HH ) 
NN  =  1 
17    CK1=HH*FCT(XX,YY) 

CK2=HH*FCT(XX+HH/2.0,YY+CK 1/2.0) 
CK3=HH*FCT(XX+HH/2.0,YY+CK2/2.0) 
CK4=HH*FCT(XX+HH,YY+CK3) 
GO  TO  (101, 102, 103), NN 

101  Yl=YY+(CKl+2.0*CK2+2.0*CK3+CK4)/6.0 
YY  =  Y1 

77    XX=XX+HH 

YEXACT=EXACT(XX) 

ERROR=YEXACT-YY 

WRITE (6, 1004)  XX,  YY,Y EX ACT, ERROR 

NN=NN+1 

IF(NN-KK)     17,17,99 

102  Y2=YY+(CKl+2.0*CK2+2.0*CK3+CK4)/6.0 
YY=Y2 

GO    TO    77 

103  Y3=YY+ (CK1+2.0~CK2+2.0«CK3+CK^ )/6.0 
YY=Y3 

GO  TO  77 
1004  FORMAT  CO'  ,9X,F15.  7  ,2X  ,F15  .7  ,16X  ,F15  .7  ,2X  ,  F15.7) 
99    RETURN 

END 

FUNCTION  FCT(XN,YN) 

FCT=-YN 

RETURN 

END 

FUNCTION  EXACTUK) 

EXACT=EXP(-XK) 

RETURN 

END 
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♦INPUT  DATA  SET  USED* 
H     =  0.5000000 
ICON  =  2 

♦♦RESULTS  OF  MILNE  PREDICTOR-CORRECTOR  METHOD** 


0.0000000 
0.5000000 
1.0000000 
1.5000000 
2.0000000 
3.0000000 
4. 0000000 
5.0000000 
6.0000000 
7.0000000 
8.0000000 
9.0000000 
10.0000000 


YP 

0.6067709 

0.3681710 

0.2233955 

0.1385594 

0.0509259 

0.0183159 

0.0061856 

0.0015049 

-0.0005340 

-0.0017356 

-0.0028179 

-0.0041233 


YC 


0.1353100 

0.0496315 

0.0181091 

0.0064749 

0.0021262 

0.0004233 

-0.0003513 

-0.0008470 

-0.0013280 


YEXACT 
1.0000000 
0.6065306 
0.367894 
0.2231301 
0.1353353 
0.0497871 
0.0183156 
0.0067379 
0.0024788 
0.0009119 
0.0003355 
0.0001234 
0.0000454 


ERROR 

-0.0002403 
-0.0002916 
-0.0002654 
0.0000253 
0.0001555 
0.0002066 
0.0002630 
0.0003525 
0.0004886 
0.0006868 
0.0009704 
0.0013734 


*INPUT  DATA  SET  USED* 
H     =  1.0000000 
ICON  =  1 

♦♦RESULTS  OF  MILNE  PREDICTOR-CORRECTOR  METHOD** 


0.0000000 
1.0000000 
2.0000000 
3.0000000 
4.0000000 
5.0000000 
6.0000000 
7.0000000 
8.0000000 
9.0000000 
10.0000000 


YP 

0.3750000 

0.1406250 

0.0527344 

0.0468752 

0.0147570 

0.0102883 

0.0014462 

-0.0005506 

-0.0039288 

-0.0064989 


YC 


0.0164931 
0.0051922 
0.0002441 
0.0005433 

-0.0009222 
0.0010040 

-0.0017257 


YEXACT 
1.0000000 
0.3678794 
0.1353353 
0.0497871 
0.0183156 
0.0067379 
0.0024788 
0.0009119 
0.0003355 
0.0001234 
0.0000454 


ERROR 

-0.0071206 

-0.0052897 

-0.0029473 

0.0018226 

0.0015457 

0.0022346 

0.0003686 

0.0012576 

-0.0008806 

0.0017711 
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PROGRAM    THREE 

NYSTROM    PREDICTOR-CORRECTOR    METHOD.    ODE:    Y»=-Y 
WIYH    Y(0)=1.    TO    SOLVE    ANOTHER    ODE    SIMPLY    CHANGE 
INITIAL    CONDITIONS    AND    FUNCTION    SUBROUTINES. 


^^  J^  ^p  *|^  »,»  *y   *^»  gfi   -»j«.  •■p.  Jf^   *^  *,»  *,*.  »j<  *f+   Jp.  -f^  *^  /,*  *■,»  *»-•  *i(*  *t*  *p.  ^j-*  *|*«  *f*  *j»  "7*  *f*   «t^  *i»  ^  *^-  «^*  ^,->.  *,»  *■,»  *»|^  ^»  >,«.  ^f»  «y»  *,-  *j-«  *y"  •<*  *l*  *P  "(*  *T"  "P  T* 

C  *    FORTRAN    PROGRAM    FOR    NYSTROM    PREDICTOR-CORRECTOR  * 

C  *    METHOD    IN    THE    NUMERICAL    SOLUTION    OF    THE    ORDINARY  * 

C  *    DIFFERENTIAL     EQUATION     IN    THE    FORM    OF    DY/DX=F(X,Y)  * 

C  *    WITH    THE     INITIAL    CONDITION    Y(X0)=Y0  * 

C  *    THE    PARAMETERS    TO    THE    PROGRAM    HAVE    THE    FOLLOWING  * 

C  *    MEANINGS...  * 

C  *         XMAX THE    TERMINAL    BOUNDARY    CONDITION  * 

C  *         EPS_ THE    CONVERGENCE    TEST    CONSTANT  * 

C  *          MAX THE    MAXIMUM    NUMBER    OF    ITERATIONS  * 

C  *         YP THE    PREDICTED    VALUE    OF    THE    NEXT  * 

C  *                                     SOLUTION    POINT  * 

C  *         YC_: THE    CORRECTED    VALUE    OF    THE    NEXT  * 

C  *                                 "SOLUTION    POINT  * 

C  *          FP THE    FUNCTION    EVALUATION  * 

C  *         H THE    STEP    SIZE    TO    BE    USED    TO    ADVANCE  * 

C  *                                    THE    POINT    OF     INTEGRATION  * 

C  *          ICON THE    DESIRED    POINT    OF    PRINTING    THE  * 

C  *                                    COMPUTED    SOLUTION    FOP     EVERY    FIXED  * 

C  *                                     NUMBER    OF    INTEGRATION    POINTS  * 

C  *         YEXACT THE    TRUE    SOLUTION    VALUES  * 

C  *         ERROR THE    ABSOLUTE    DEVIATION    OF    THE    NUMERICAL* 

C  *                                      SOLUTION    FROM    THE    TRUE    SOLUTION  * 

C  *         FCT FUNCTION    SUBROUTINE    TO    EVALUATE    ODE  * 

C  *                                    FUNCTION.  * 

C  *         EXACT FUNCTION    SUBROUTINE    TO    COMPUTE    TRUE  * 

C  *                                  "SOLUTION    VALUES  * 

C  *         RKUTTA THE    SUBROUTINE    USED    TO    GENERATE  * 

C  *                                    NEEDED    STARTING    VALUES  * 

£  ****************************************************** 

C 
C 

C  INITIALIZE    INPUT    CONSTANTS 
C 

XMAX=10.0 

MAX  =  2 

EPS=0. OOOOOl 

c 

C  READ    STEP    SIZE    AND    CONTROL    PRINTING 


5  READ{5,100)    H,ICON 

C 
C  SPECIFY    INITIAL    CONDITION 

X0=0.0 

Y0=1.0 

YR=Y0 
C 

C  TEST    FOR    END    OF     INPUT    DATA 

C 

IF(ICON-O)    45,45,7 
C 

C  PRINT    INPUT    DATA    SET 

C 
7  WRITE(6,1000)    H, ICON 

WRITE{6,1001)     X0,Y0 
C 

C  COMPUTE    NEXT    SOLUTION    POINT 

C 
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NC=( XMAX+H/2.)/H 

X1=X0+H 

X=X0+2.*H 

CALL    RKUTTA(XO,YR, Y1,H) 

F1=FCT(X1,Y1) 

DO    40    N=2,NC 

M=l 

YP=Y0+2.*H*F1 

FP=FCT(X,YP) 
11         YC=Y1+H/2.*(F1+FP) 

DELY  =  ABS(YOYP) 
C 
C  TEST    FOR    CONVERGENCE 


IF(DELY-EPS)    30,30^,15 
C 
C  TEST     IF    MAXIMUM    NUMBER    OF    ITERATIONS    IS    SATISFIED 


15    IF(M-MAX)  20,20,30 
20    FP=FCT(X,YC) 
M=M+1 
GO  TO  11 
C 

C      COMPUTE  TRUE  SOLUTION  VALUE 
C 

30    YEXACT=EXACT(X) 
C 

C      COMPUTE  ERROR  IN  CALCULATION 
C 

ERROR=YEXACT-YC 
C 

C      TEST  IF  DESIRED  POINT  OF  PRINTING  IS  REACHED 
C 

IF( (N/ICON*ICON).EQ.N)  WRIT E (6 , 1 002 )  X , YP, YC , YEXACT , 
1ERR0R 
C 

C      UPDATE  REQUIRED  PARAMETER  VALUES 
C 

Y0=Y1 
Y1=YC 
F1=FP 
X  =  X  +  H 
40  CONTINUE 
C 

C      LOOP  BACK  TO  READ  NEXT  SET  OF  INPUT  DATA 
C 

GO  TO  5 
45    STOP 
100   FORMAT(F10*7,I10) 

1000  F0RMAT{////41X,f*INPUT  DATA  SET  US ED*» //43X , «H     =', 
1F10.7/43X, •ICON  =•  ,  I  10) 

1001  FORMAT (//29X»,*RESULTS  OF  NYSTROM  PREDICTOR-CORRECTOR 
1METH0D  *,//10Xf«    X   «,5X,'    YP    ',5X,'    YC    ', 

15X,«    YEXACT    ',5X,»    ERROR    • //12X, F10. 7 ,44X , 
1F10.7) 

1002  FORMAT ( '0»,9X,F15.7,2X*F15.7,2X,F15.7,2X,F15.7,2X, 
1F15.7) 

END 

SUBROUTINE  RKUTTA( XX ,YY, Yl , HH ) 
CK1=HH*FCT(XX,YY) 

CK2=HH*FCT(XX+HH/2.0,YY+CK 1/2.0) 
CK3=HH*FCT(XX+HH/2.0,YY+CK2/2iO) 
CK4=HH*FCT(XX+HH,YY+CK3) 
Yl=YY+(CKl+2.0*CK2+2.0*CK3+CK4)/6.0 
XX=XX+HH 
YEXACT=EXACT(XX) 
ERR0R=YEXACT-Y1 

WRITE (6, 1004)  XX, Yl, YEXACT, ERROR 
1004  FORMAT (  ,0«,9X,F15.7,2X,F15.7,16X,F15.7,2X,F15.7) 
RETURN 
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END 


FUNCTION 
FCT=-YN 
RETURN 
END 


FCT(XN,YN) 


FUNCTION  EXACT(XK) 

EXACT=EXP(-XK) 

RETURN 

END 
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*INPUT  DATA  SET  USED* 
H     =  0.5000000 
ICON  =  2 

^RESULTS  OF  NYSTROM  PREDICTOR-CORRECTOR  METHOD  * 


0.0000000 
0.5000000 
1.0000000 
2.0000000 
3.0000000 
4.0000000 
5.0000000 
6.0000000 
7.0000000 
8.0000000 
9.0000000 
10.0000000 


YP 

0.60677  09 
0.3932291 
0.1444499 
0.0516075 
0.0184066 
0.0065643 
0.0023410 
0.0008349 
0.0002977 
0.0001062 
0.0000379 


YC 


0.3636068 
0.1298206 
0.0463004 
0.0165119 
0.0058886 
0.0021000 
0.0007489 
0.0002671 
0.0000952 
0.0000340 


YEXACT 
1.0000000 
0.6065306 
0.3678794 
0.1353353 
0.0497871 
0.0183156 
0.0067379 
0.0024788 
0.0009119 
0.0003355 
0.0001234 
0.0000454 


ERROR 

-0.0002403 
0.0042726 
0.0055147 
0.0034867 
0.0018037 
0.0008494 
0.0003788 
0.0001630 
0.0000684 
0.0000282 
0.0000114 


*INPUT    DATA    SET    USED* 
H  =    1.0000000 

ICON    =  1 

^RESULTS    OF    NYSTROM    PREDICTOR-CORRECTOR    METHOD    * 


0.0000000 
1.0000000 
2.0000000 
3.0000000 
4.0000000 
5.0000000 
6.0000000 
7.0000000 
8.0000000 
9.0000000 
10.0000000 


YP 

0.3750000 
0.2500000 
0. 0625000 
0.0468750 

-0.0078125 
0.0097656 

-0.0087891 
0.0046387 

-0.0045166 
0.0028381 


YC 


0.1093750 
0.0156250 
-0.0053594 
-0.0078125 
-0.0041504 
-0.0021973 
-0.0005798 
-0.0003052 
0.0000572 


YEXACT 
1.0000000 
0.3678794 
0.1353353 
0.0497871 
0.0183156 
0.0067379 
0.0024788 
0.0009119 
0.0003355 
0.0001234 
0.0000454 


ERROR 

-0.0071206 
0.0259603 
0.0341621 
0.0241750 
0.0145504 
0.0066291 
0.0031091 
0.0009153 
0.0004286 

-0.0000118 
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PROGRAM    FOUR 

HERMITE    PREDICTOR-CORRECTOR    METHOD.    ODE:    Y«=-Y 
WIYH    Y(0)=1.    TO    SOLVE    ANOTHER    ODE    SIMPLY    CHANGE 
INITIAL    CONDITIONS    AND    FUNCTION    SUBROUTINES. 


C*V  <JU  *JU  <jU  -^  Vf  J    .'.  f  Jv.'/  J.  vi.  Vr  V'  *''  ^f  *)*  V?  -»u  -1'  l'-  i'*  -1*  °  -  -■**  ***  -1'  J,  J'  4,  »'.  J.  J.-  J-  j,  o,  j,  j-  o-  J/  u,  w.  O.  J.  ^  *),  g-  v^  »i.  JU  J.  ^  *i. 
*^  *?  "P  *r  *»*  nr  *.*  ^r  *r-  *r*  *<*  v  *"»*  *r  3<*  *r*  *tt  o*  ■nr  "*  •*?•  *n  nr*  *ir  »■"*  i-  '<*  "V*  *v*  ***  *Tr  t*  *r- *v  **»  ttttt^T'C*  v  v^«  ^tt'pt  *"*■  *r  n* 

C  *    FORTRAN    PROGRAM    FOR    HERMITE    PREDICTOR-CORRECTOR  * 

C  *    METHOD    IN    THE    NUMERICAL    SOLUTION    OF    THE    ORDINARY  * 

C  *    DIFFERENTIAL    EQUATION    IN    THE    FORM    OF    DY/DX=F(X,Y)  * 

C  *    WITH    THE    INITIAL    CONDITION    Y(XO)=YO  * 

C  *    THE    PARAMETERS    TO    THE    PROGRAM    HAVE    THE    FOLLOWING  * 

C  *    MEANINGS...  * 

C  *         XMAX THE    TERMINAL    BOUNDARY    CONDITION  * 

C  *         EPS THE    CONVERGENCE    TEST    CONSTANT  * 

C  *         MAX , THE    MAXIMUM    NUMBER    OF    ITERATIONS  * 

C  *         YP THE    PREDICTED    VALUE    OF    THE    NEXT  * 

C  *  SOLUTION    POINT  * 

C  *         YC THE    CORRECTED    VALUE    OF    THE    NEXT  * 

C  *  SOLUTION    POINT  * 

C  *         FP THE    FUNCTION    EVALUATION  * 

C  *         H THE    STEP    SIZE    TO    BE    USED    TO    ADVANCE  * 

C  *  THE    POINT    OF     INTEGRATION  * 

C  *  ICON THE    DESIRED    POINT    OF    PRINTING    THE  * 

C  *  COMPUTED    SOLUTION    FOR    EVERY    FIXED  * 

C  *  NUMBER    OF     INTEGRATION    POINTS  * 

C  *         YEXACT THE    TRUE    SOLUTION    VALUES  * 

C  *         ERROR THE    ABSOLUTE    DEVIATION    OF    THE    NUMERICAL* 

C  *  SOLUTION    FROM    THE    TRUE    SOLUTION  * 

C  *         FCT FUNCTION    SUBROUTINE    TO    EVALUATE    ODE 

C  *  FUNCTION  * 

C  *         EXACT FUNCTION    SUBROUTINE    TO    COMPUTE    TRUE  * 

C  *  SOLUTION    VALUES  * 

C  *  RKUTTA THE    SUEROUTINE    USED    TO    GENERATE  * 

C  *  NEEDED    STARTING    VALUES  * 

C 
C 

C  INITIALIZE    INPUT    CONSTANTS 
C 

XMAX=10.0 

MAX=2 

EPS=0. 000001 
C 

C  READ    STEP    SIZE    AND    CONTROL    PRINTING 


C 

i 

C 

c 

c 


5  READ(5,100)    H,IC0N 

C 

SPECIFY    INITIAL    CONDITION 

X0=0.0 

Y0=1.0 

YR=Y0 
C 

C  TEST    FOR    END    OF    INPUT    DATA 

C 

IF(ICON-O)    45,45,7 
C 

C  PRINT     INPUT    DATA    SET 

C 
7  WRITE(6,1000)    H,ICON 

WRITE16, 1001 )    X0,Y0 
C 

C  COMPUTE    NEXT    SOLUTION    POINT 

C 
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NC=(XMAX+H/2.)/H 

X1=X0+H 

X=X0+2.*H 

CALL  RKUTTA(X0,YR,Y1,H) 

F0=FCTIX0,Y0) 

F1=FCT{X1,Y1) 

DO  40  N=2,NC 

M=l 

YP=5.*Y0-4.*Y1+H*(4.*F1+2.*F0) 

FP=FCT(X,YP) 
11    YC=Y0+H/3.*(F0+4.*F1+FP) 

DELY=ABS(YC-YP) 
C 

C      TEST  FOR  CONVERGENCE 
C 

IF(DELY-EPS)  30,30,15 
C 

C      TEST  IF  MAXIMUM  NUMBER  OF  ITERATIONS  IS  SATISFIED 
C 
15    IF(M-MAX)  20,20,30 
20    FP=FCT(X,YC) 

M=M+1 

GO  TO  11 
C 

C      COMPUTE  TRUE  SOLUTION  VALUE 
C 

30    YEXACT=EXACT(X) 
C 

C      COMPUTE  ERROR  IN  CALCULATION 
C 

ERROR=YEXACT-YC 
C 

C  TEST    IF    DESIRED    POINT    CF    PRINTING    IS    REACHED 

C 

IF( (N/ICON*ICON).EQ.N)    WRITE (6, 1002 )    X, YP, YC , YEXACT, 
1ERR0R 
C 

C  UPDATE    REQUIRED    PARAMETER    VALUES 

C 

Y0=Y1 

Y1=YC 

F0=F1 

F1  =  FP 

X=X  +  H 
40         CONTINUE 
C 

C  LOOP    BACK    TO    READ    NEXT    SET    OF    INPUT    DATA 

C 

GO    TO    5 
45  STOP 

100       F0RMAT(F10.7,I10) 

1000  F0RMAT(////41X,»*INPUT    DATA    SET    USED* • //43X , • H  =', 
1F1C.7/43X,1 ICON    =«,I10) 

1001  F0RMAT(//29X,»*RESULTS    OF    HERMITE    PREDICTOR^CORP ECTOR 
1METH0D    *«//10Xf«  X       • ,5X,«  YP  • ,5X,«  YC  ', 

15X,'  YEXACT  ,,5X,'  ERROR  • //12X, F 10. 7 ,44X, 

1F10.7) 

1002  FORMAT (  •0" »9Xf F15. 7t 2XtF15.7f 2X,F15. 7f 2X,F15.7,2X , 
1F15.7) 

END 

SUBROUTINE  RKUTT A( XX , YY , Yl , HH ) 
CK1=HH*FCT(XX,YY) 

CK2=HH*FCT(XX+HH/2.0,YY+CK 1/2.0) 
CK3=HH*FCT(XX*HH/2.0,YY+CK2/2.0J 

CK4=HH*FCT(XX+HH,YY+CK3) 

Yl=YY+(CKl+2.0*CK2+2.0*CK3+CK4)/6.0 

XX=XX+HH 

YEXACT=EXACT(XX) 

ERR0R=YEXACT~Y1 

WRITE(6,1004)  XX, Yl, YEXACT, ERROR 
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1004  FORMAT ( «0',9X,F15.7,2X,F15.7,16X,F15.7,2XfF15.7) 
RETURN 
END 

FUNCTION  FCT(XN,YN) 

FCT=-YN 

RETURN 

END 

FUNCTION  EXACT(XK) 

EXACT=EXPt-XK) 

RETURN 

END 
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*INPUT  DATA  SET  USED* 
H     =  0.5000000 
ICON  =  2 

^RESULTS  OF  HERMITE  PREDICTOR-CORRECTOR  METHOD  * 


0.0000000 
0.5000000 
1.0000000 
2.0000000 
3.0000000 
4.0000000 
5.0000000 
6.0000000 
7.0000000 
8.0000000 
9.0000000 
10.0000000 


YP 

0.6067709 

0.3593760 

0.1296880 

0.0439692 

0.0106042 

-0.0043826 

-0.0139239 

-0.0234190 

-0.0358057 

-0.0535792 

-0.0797584 


YC 


0.3675976 

0.1349390 

0.0491728 

0.0173813 

0.0053371 

0.0003908 

-0.0021938 

-0.0042811 

-0.0067376 

-0.0101509 


YEXACT 
1.0000000 
0.6065306 
0.3678794 
0.1353353 
0.0497871 
0.0183156 
0.0067379 
0.0024788 
0.0009119 
0.0003355 
0.0001234 
0.0000454 


ERROR 

-0.0002403 
0.0002818 
0.0003963 
0.0006143 
0.0009343 
0.0014008 
0.0020879 
0.0031057 
0.0046165 
0.0068610 
0.0101963 


*INPUT  DATA  SET  USED* 
H     =  1.0000000 
ICON  =  1 

^RESULTS  OF  HERMITE  PREDICTOR-CORRECTOR  METHOD  * 


0.0000000 
1.0000000 
2.0000000 
3.0000000 
4.0000000 
5.0000000 
6.0000000 
7.0000000 
8.0000000 
9.0000000 
10.0000000 


YP 

0.3750000 
0.0000000 
0.1620350 

-0.2105594 
0.3834665 

0.6344536 
1.0731880 

-1.8065000 
3.0441850 

-5.1285590 


YC 


0.1296299 

0.0732166 
-0.0092715 

0. 0599074 
-0.0839149 

0.1479157 
-0.2466852 

0.4165850 
-0.7014816 


YEXACT 
1.0000000 
0.3678794 
0.1353353 
0.0497871 
0.0183156 
0.0067379 
0.0024788 
0.0009119 
0.0003355 
0.0001234 
0.0000454 


ERROR 

-0.0071206 

0.0057054 
-0.0234295 

0.0275871 
-0.0531695 

0.0863936 
-0.1470038 

0.2470207 
-0.4164615 

0.7015270 
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PROGRAM    FIVE 

HAMMING    PREDICTOR-CORRECTOR    METHOD.    ODE:     Y"=-Y 
WIYH    Y(0)=1.    TO    SOLVE    ANOTHER    ODE    SIMPLY    CHANGE 
INITIAL    CONDITIONS    AND    FUNCTION    SUBROUTINES. 


■*(-*  rg*  wfr  -"<^  ^^  '.*  -*»*■  *»-  "j~  -*»"  -i*  ?.K  "i*  *j*  1*  **»*  *$•*?•  ^r  ^^^^^vov^T*^^*^^*iv'r'r''i*'r'ii*r'i'''t*'T,'T'*i'¥'r'rn**r'  ~&  *fr  *r-  t*  *c- 

C              *    FORTRAN    PROGRAM    FOR    HAMMING    PREDICTOR-CORRECTOR  * 

C               *    METHOD    IN    THE    NUMERICAL    SOLUTION    OF    THE    ORDINARY  * 

C               *    DIFFERENTIAL    EQUATION     IN    THE    FORM    OF    DY/DX=F(XtY)  * 

C               *    WITH    THE    INITIAL    CONDITION    Y(X0)=Y0  * 

C               *    THE    PARAMETERS    TO    THE    PROGRAM    HAVE    THE    FOLLOWING  * 

C               *    MEANINGS...  * 

C               *         XMAX THE    TERMINAL    BOUNDARY    CONDITION  * 

C               *         EPS THE    CONVERGENCE    TEST    CONSTANT  * 

C               *         MAX THE    MAXIMUM    NUMBER    OF    ITERATIONS  * 

C               *         YPM Z PREDICTED    SOLUTION    VALUE  * 

C               *          YP    MODIFIED    PREDICTED    SOLUTION    VALUE  * 

C               *         YCffi CORRECTED    SOLUTION    VALUE  * 

C               *         YC .^MODIFIED    CORRECTED    SOLUTION    VALUE  * 

C               *          FP THE    FUNCTION    EVALUATION  * 

C               *         H THE    STEP    SIZE    TO    BE    USED    TO    ADVANCE  * 

C               *                                    THE    POINT    OF     INTEGRATION  * 

C               *          ICON THE    DESIRED    POINT    OF    PRINTING    THE  * 

C                *                                      COMPUTED    SOLUTION    FOR    EVERY     FIXED  * 

C               *                                    NUMBER    OF    INTEGRATION    POINTS  * 

C               *         YEXACT THE    TRUE    SOLUTION    VALUES  * 

C  *         ERROR THE    ABSOLUTE    DEVIATION    OF    THE    NUMERICAL* 

C               *                                      SOLUTION    FROM    THE     TRUE    SOLUTION  * 

C               *         FCT FUNCTION    SUBROUTINE    TO    EVALUATE    ODE  * 

C               *                                    FUNCTION  * 

C               *         EXACT FUNCTION    SUBROUTINE    TO    COMPUTE    TRUE  * 

C               *                                     SOLUTION    VALUES  * 

C               *         RKUTTA THE    SUBROUTINE    USED    TO    GENERATE  * 

C               *                                    NEEDED    STARTING    VALUES  * 

C5^  -i  -V  -V  *V  ^  ^-  f  ''•  l'-  '••  ^  -,J  *''  t1-  -1'  ^  ^-  ^  *'•  v  •'  -'--'-  U-  .'--■-,'.  v-  *J-    -■-  o.  -J-  >  -  -  -  0 .  .~  -"-  o.  y-  ^y  v.  ,*,  . V  >^  V-  ,v,  V*  V-  ~«-  J»  ■«*-  -*-  -<- 
t*  *r  **h  *r   Jr  *v  *r-  *f  fp  "V-  -v  *r*  *v  «v»  •*&  f-  n*  *r*  *r*  •¥■  nr  *<-  i*  "i**  *i-  ^r  *»-***  *t*  n*  "i*  t*  ***  **■*  "*»"•  "V  *r*  *v»  *i*  t*  ^i*1  "i*  ***■  ^i*  n*  <t*  *nr  *t*  n*  -r>  *■**  -»-  */•  *f* 

C 
C 

C  INITIALIZE    INPUT    CONSTANTS 

C 

XMAX=10.0 

MAX-2 

EPS=0. 000001 

K=3 
C 
C  READ    STEP    SIZE    AND    CONTROL    PRINTING 


C 
i 

C' 

c 
c 


5  READ(5,100)    H,ICON 

C 

SPECIFY    INITIAL    CONDITION 

ET-0.0 

X0=0.0 

Y0=1  .0 

YR  =  YO 
C 

C  TEST    FOR    END    OF    INPUT    DATA 

C 

IF(ICON-O)    45,45,7 
C 
C  PRINT     INPUT    DATA    SET 


C 


7  WRITE(6, 1000)    H,ICON 

WRITE(6,1001 )    XO,YO 
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C      COMPUTE  NEXT  SOLUTION  POINT 
C 

NC=(XMAX+H/2.)/H 

X1=X0+H 

X2=X0+2.*H 

X3=X0+3.*H 

X=X0+4.*H 

CALL  RKUTTA(X0,YR,Y1,Y2,Y3,K,H) 

F1=FCT(X1,Y1) 

F2=FCT(X2,Y2) 

F3=FCT(X3tY3) 

DO  40  N=4,NC 

M=l 

YPM=Y0+4.*H/3.*(2.*F1-F2+2.*F3) 

YP=YPM-112./121.*ET 

FP=FCT(X,YP) 
11    YCM=1./8.*(9.*Y3-Y1 ) +3 .*H/ 8.*< FP+2. *F3-F2) 

ET=YPM-YCM 

YC=YCM+9./121.*ET 

DELY=ABS(YC-YP) 
C 

C      TEST  FOR  CONVERGENCE 
C 

IF(DELY-EPS)  30,30,15 
C 

C  TEST    IF    MAXIMUM    NUMBER    OF    ITERATIONS    IS    SATISFIED 

C 
15  IF(M-MAX)     20,20,30 

20         FP=FCT(X,YC) 

M=M+1 

GO    TO    11 
C 

C      COMPUTE  TRUE  SOLUTION  VALUE 
C 

30    YEXACT=EXACT{X) 
C 

C      COMPUTE  ERROR  IN  CALCULATION 
C 

ERROR=YEXACT-YC 
C 

C      TEST  IF  DESIRED  POINT  OF  PRINTING  IS  REACHED 
C 

IF( (N/ICON*ICON).EQ.N)  WRITE (6 , 1002 )  X , YP, YC , YEXACT , 
1ERR0R 
C 

C      UPDATE  REQUIRED  PARAMETER  VALUES 
C 

F1  =  F2 

F2  =  F3 

F3=FP 

Y0=Y1 

Y1=Y2 

Y2=Y3 

Y3=YC 

X  =  X  +  H 
40    CONTINUE 
C 

C      LOOP  BACK  TO  READ  NEXT  SET  OF  INPUT  DATA 
C 

GO    TO    5 
45         STOP 
100       F0RMAT(F10.7,I10) 

1000  F0RMATC////41X, »*INPUT    DATA    SET    USED*1 //43X , • H  =■ , 
1F10.7/43X, ' ICON    =*,I10) 

1001  F0RMAT(//29X,'*RESULTS    OF    HAMMING    PREDICTOR-CORRECTOR 
1METH0D    *«//10X,«  X       «,5X,«  YP  ',5X,«  YC  ', 

15X,»  YEXACT  '^X,'  ERROR  • //12X,  F10.7 ,44X, 

1F10.7) 

1002  FORMAT i  »0» f9X, F15.7f 2XtF15,7f 2Xf F15.7, 2X,F15.7,2X, 
1F15.7) 

END 
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SUBROUTINE  RKUTT A( XX , YY, Yl , Y2, Y3 »KK, HH ) 

NN=1 
17    CK1=HH*FCT(XX,YY) 

CK2=HH*FCT(XX+HH/2.0,YY+CKl/2.0) 

CK3=HH*FCT(XX+HH/2.0,YY+CK2/2#0) 

CK4=HH*FCT(XX+HH,YY+CK3) 

GO  TO  (101, 102, 103), NN 
101   Yl=YY+(CKl+2.0*CK2+2.0*CK3+CK4)/6.0 

YY=Y1 
77    XX=XX+HH 

YEXACT=EXACT(XX) 

ERROR=YEXACT-YY 

WRITE(6,1004)  XX, YY,YEXACT, ERROR 

NN=NN+1 

IF(NN-KK)  17,17,99 
10  2   Y2=YY+(CKl+2.0*CK2+2.0*CK3+CK4)/6. 0 

YY  =  Y2 

GO  TO  77 
103   Y3=YY+(CKl+2.0*CK2+2.0*CK3+CK4)/6. 0 

YY=Y3 

GO  TO  77 
1004  FORMAT ( «0» ,9X, Fl 5. 7, 2X ,F 15. 7 ,16X ,F 15.7 ,2X, F15.7) 
99    RETURN 

END 

FUNCTION  FCT(XN,YN) 

FCT=-YN 

RETURN 

END 

FUNCTION  EXACT(XK) 
EXACT=EXP(-XK) 

RETURN 
END 
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♦INPUT  DATA  SET  USED* 
H     =  0.5000000 
ICON  =  2 

♦RESULTS  OF  HAMMING  PREDICTOR-CORRECTOR  METHOD  * 


0.0000000 
0.5000000 
1.0000000 
1.5000000 
2.0000000 
3.0000000 
4.0000000 
5.0000000 
6.0000000 
7.0000000 
8.0000000 
9.0000000 
10.0000000 


YP 

0.6067709 
0.3681710 
0.2233955 
0.1385594 
0.0494538 
0.0184383 
0.0068499 
0.0025701 
0.0009748 
0.0003766 
0.0001496 
0.0000619 


YC 


0.1355409 
0.0499289 
0.0183963 
0.0067794 
0.0024990 
0.0009215 
0.0003401 
0.0001256 
0.0000465 


YEXACT 
1.0000000 
0.6065306 
0.3678794 
0.2231301 
0.1353353 
0.0497871 
0.0183156 
0.0067379 
0.0024788 
0.0009119 
0.0003355 
0.0001234 
0.0000454 


ERROR 

-0.0002403 
-0.0002916 
-0.0002654 
-0.0002056 
-0.0001418 
-0.0000806 
-0.0000415 
-0.0000202 
-0.0000096 
-0.0000046 
-0.0000022 
-0.0000011 


♦INPUT  DATA  SET  USED* 
H     =  1.0000000 
ICON  =  1 

♦RESULTS  OF  HAMMING  PREDICTOR-CORRECTOR  METHOD  * 


0.0000000 
1.0000000 
2.0000000 
3.0000000 
4.0000000 
5.0000000 
6.0000000 
7.0000000 
8.0000000 
9.0000000 
10.0000000 


YP 

0.3750000 
0.1406250 
0.0527344 
0.0468752 

-0.0199183 
0.0245465 

-0.0516121 
0.0793315 

-0.1184798 
0.1838089 


YC 


0.0190868 
0.0056581 
0.0057380 

-0.0008928 
0.0053606 

-0.0069360 
0.0110968 


YEXACT 
1.0000000 
0.3678794 
0.1353353 
0.0497871 
0.0183156 
0.0067379 
0.0024788 
0.0009119 
0.0003355 
0.0001234 
0.0000454 


ERROR 

-0.0071206 
-0.0052897 
-0.0029473 
-0.0007712 

0.0010798 
-0.0032592 

0.0018047 
-0.0050251 

0.0070594 
-0.0110514 
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c 

I 

c' 
c 
c 


PROGRAM    SIX 

SECOND    ORDER    ADAMS    PREDICTOR-CORRECTOR    METHOD.    ODE 
Y«=-Y    WITH    Y(0)=1.    TO    SOLVE    ANOTHER    ODE    CHANGE 
INITIAL    CONDITIONS    AND    FUNCTION    SUBROUTINES. 


Cy,  J-  J-  u*  y,  .1,  yL,  -.c  -t-  ■*!.  »j,  j,  of  -J*  *>.  o-  y-  -  ,«-  v  -  o.  -jl.  .j-  j-  -J,  o,  a-  J  -  J^  a-  ol-  -j-  J-  -j-  J-  y-  «*"  Vt-  "*-  V-  *■»-  **-  >*~  \V  *■*-  V-  V-  ••*-  -*•  -J*  *>**  V~  -'-  -'-  -'' 
•P  i*  *p  *i*  *r-  *v  t*  t*  *»■*  *>"  *r*  *v  n*  *"r  "<*  ***  n*  t*  hp  *v  ■*?*  "r-  *f  t*  *v  "nr  ■v*  *t*  *r*  *f  n**  ^h  *f*  ■?(*  *r-  *r*  *r-  *f  *»*  t*  n*  *r*  ^r  *>*  Jr  *r  *r-  *r  -i*  t*  *r*  *»*  "V*  *i* 

C  *    FORTRAN    PROGRAM    FOR    SECOND    ORDER    ADAMS    PREDICTOR-  * 

C  *    CORRECTOR    METHOD    IN    THE    NUMERICAL    SOLUTION    OF    THE  * 

C  *    ORDINARY    DIFFERENTIAL    EQUATION    IN    THE    FORM    OF  * 

C  *    DY/DX=F(X,Y)     WITH    THE    INITIAL    CONDITION    Y(X0)=Y0  * 

C  *    THE    PARAMETERS    TO    THE    PROGRAM    HAVE    THE    FOLLOWING  * 

C  *    MEANINGS...  * 

C  *         XMAX THE    TERMINAL    BOUNDARY    CONDITION  * 

C  *          EPS_ THE    CONVERGENCE    TEST    CONSTANT  * 

C  *         MAX____ THE    MAXIMUM    NUMBER    OF    ITERATIONS  * 

C  *         YP THE    PREDICTED    VALUE    OF    THE    NEXT  * 

C  *                                     SOLUTION    POINT  * 

C  *         YC THE    CORRECTED    VALUE    OF    THE    NEXT  * 

C  *                          "^        SOLUTION    POINT  * 

C  *         FP THE    FUNCTION    EVALUATION  * 

C  *         H THE    STEP    SIZE    TO    BE    USED    TO    ADVANCE  * 

C  *                                    THE    POINT    OF     INTEGRATION  * 

C  *          ICCN_              THE    DESIRED    POINT    OF    PRINTING    THE  * 

C  *                                    COMPUTED    SOLUTION    FOP    EVERY    FIXEO  * 

C  *                                     NUMBER    OF    INTEGRATION    POINTS  * 

C  *         YEXACT THE    TRUE    SOLUTION    VALUES  * 

C  *  ERROR THE    ABSOLUTE    DEVIATION    OF    THE    NUMERICAL* 

C  *                                     SOLUTION    FROM    THE    TRUE    SOLUTION  * 

C  *         FCT FUNCTION    SUBROUTINE    TO    EVALUATE    ODE  * 

C  *                                    FUNCTION.  * 

C  *         EXACT FUNCTION    SUBROUTINE    TO    COMPUTE    TRUE  * 

C  *                                     SOLUTION    VALUES  * 

C  *         RKUTTA THE    SUBROUTINE    USED    TO    GENERATE  * 

C  *                                     NEEDED    STARTING    VALUES  * 

C 

C 

C  INITIALIZE    INPUT    CONSTANTS 
C 

XMAX=10.0 

MAX=2 

EPS=0. OOOOOl 

c 

C  READ    STEP    SIZE    AND    CONTROL    PRINTING 


5  READ(5,100)    H, ICON 

C 

SPECIFY    INITIAL    CONDITION 

X0=0.0 
Y 0*1.0 

YR=YO 
C 

C  TEST    FOR    END    OF    INPUT    DATA 

C 

IF(ICON-O)    45,45,7 
C 

C  PRINT     INPUT    DATA    SET 

C 
7  WRITE( 6,1000)    H,ICON 

WRITE16, 1001)    X0,Y0 
C 

C  COMPUTE    NEXT    SOLUTION    POINT 

C 
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NC=(XMAX+H/2.)/H 

FO=FCT(XO,YO) 

X1=X0+H 

X=X0+2.*H 

CALL    RKUTTA(X0,YR,Y1,H) 

F1=FCT(X1,Y1) 

DO   40    N=2,NC 

M=l 

YP=Y1+H/2.*(3.*F1-F0) 

FP=FCT(X,YP) 
11         YC=Y1+H/2.*(FP+F1) 

DELY=ABS(YC-YP) 
C 

C  TEST    FOR    CONVERGENCE 

C 

IF(DELY-EPS)    30,30,15 
C 

C      TEST  IF  MAXIMUM  NUMBER  OF  ITERATIONS  IS  SATISFIED 
C 
15    IF(M-MAX)  20,20,30 
20    FP=FCT(X,YC) 

M=M+1 

GO  TO  11 
C 

C      COMPUTE  TRUE  SOLUTION  VALUE 
C 

30    YEXACT=EXACT(X) 
C 

C      COMPUTE  ERROR  IN  CALCULATION 
C 

ERROR=YEXACT-YC 
C 

C      TEST  IF  DESIRED  POINT  CF  PRINTING  IS  REACHED 
C 

IF( (N/ICON*ICON).EQ.N)  WRITE (6 , 1002)  X, YP , YC , YEXACT, 
1ERR0R 
C 

C      UPDATE  REQUIRED  PARAMETER  VALUES 
C 

Y0=Y1 

Y1=YC 

F1=FP 

y  —  y  j.u 

40         CONTINUE 

r 

C  LOOP    BACK    TO    READ    NEXT    SET    OF    INPUT    DATA 

C 

GO    TO    5 

45         STOP 

100       FORMAT(F10.7,I10) 

1000  F0RMAT{////41X,«*INPUT    DATA    SET    US ED*« //43X ,  •  H  =  • 
1F10.7/43X,»ICON    =»,I10) 

1001  FORMAT (//29X,«**RESULTS    OF    ADAM2    PREDICTOR-CORRECTOR 
1METHCD**'//10X, ■  X       ' ,5X,»  YP  ' ,5X,»  YC  ' 

15X,'  YEXACT  ',5X,'  ERROR  • //12 X, F10.7  ,44X , 

1F10.7) 

1002  FORMAT("0«  ,9X ,F15.7 , 2X , F15 .7 ,2 X , F15 .7 , 2X , F15 .7, 2X , 
1F15.7) 

END 

SUBROUTINE  RKUTTAt XX , YY , Yl ,HH ) 
CK1=HH*FCT{XX,YY) 

CK2=HH*FCT(XX+HH/2.0,YY+CK1/2*0J 
CK3=HH*FCT(XX+HH/2.0,YY+CK2/2»0) 
CK4=HH*FCT(XX+HH,YY+CK3) 
Yl=YY+(CKl+2.0-CK2+2.0*CK3+CK4)/6.0 
XX=XX+HH 

YEXACT=EXACT(XX} 
ERR0R=YEXACT-Y1 

WRITE (6, 1004)  XX, Yl, YEXACT, ERROR 
1004  FORMAT (,0,,9X,F15.7,2X,F15.7,16X,F15.7,2X,F15.7) 
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RETURN 
END 


FUNCTION 
FCT=-YN 

RETURN 
END 


FCT(XN,YN) 


FUNCTION  EXACT(XK) 

EXACT=EXP(-XK) 

RETURN 

END 
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*INPUT    DATA    SET    USED* 
H  =    0.5000000 

ICON    =  2 

**RESULTS    OF    ADAM2    PREDICTOR-CORRECTOR    METHOD** 


0.0000000 
0.5000000 
1.0000000 
2.0000000 
3.0000000 
4.0000000 
5.0000000 
6.0000000 
7.0000000 
8.0000000 
9.0000000 
10.0000000 


YP 

0.6067709 
0.4016927 
0.2968013 
0.2556222 
0.2401217 
0.2342887 
0.2320936 
0.2312676 
0.2309567 
0.2308398 
0.2307958 


YC 


0.3634746 

0.1248231 

0.0349784 

0.0011687 

-0.0115542 

-0.0163420 

-0.0181437 

-0.0188217 

-0.0190768 

-0.0191728 


YEXACT 
1.0000000 
0.6065306 
0.3678794 
0.1353353 
0.0497871 
0.0183156 
0.0067379 
0.0024788 
0.0009119 
0.0003355 
0.0001234 
0.0000454 


ERROR 

-0.0002403 
0.0044048 
0.0105122 
0.0148087 
0.0171469 
0.0182921 
0.0188207 
0.0190556 
0.0191571 
0.0192002 
0.0192182 


*INPUT  DATA  SET  USED* 
H     =  1.0000000 
ICON  =  1 

**RESULTS  OF  ADAM2  PREDICTOR-CORRECTOR  METHOD** 


0.0000000 
1.0000000 
2.0000000 
3.0000000 
4.0000000 
5.0000000 
6.0000000 
7.0000000 
8.0000000 
9.0000000 
10.0000000 


YP 

0.3750000 
0.3125000 
0.3437500 
0.3281250 
0.3359375 
0.3320313 
0.3339844 
0.3330078 
0.3334961 
0.3332520 


YC 


0.1015625 
-0.0312500 
-0.0996094 
-0.1328125 
-rO. 1499023 
-0.1582031 
-0.1624756 
-0.1645508 
-0.1656189 


YEXACT 
1.0000000 
0.3678794 
0.1353353 
0.0497871 
0.0183156 
0.0067379 
0.0024788 
0.0009119 
0.0003355 
0.0001234 
0.0000454 


ERROR 

-0.0071206 
0.0337728 
0.0810370 
0.1179250 
0.1395504 
0.1523811 
0.1591150 
0.1628110 
0.1646742 
0.1656643 
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PROGRAM    SEVEN 

THIRD    ORDER    ADAMS    PREDICTOR-CORRECTOR    METHOD.    ODE 
Y*=-Y    WITH    Y(0)=1.    TO    SOLVE    ANOTHER    ODE    CHANGE 
INITIAL    CONDITIONS    AND    FUNCTION    SUBROUTINES. 


C«V  «**  •**  *■  V*  V-  V*  -A-  Vr  -''  °r  ^-  -J*  **'  *■**  *  *f*  V*  *J*  "^  *fe  **-•  ***  ***  ■**■*  A  **■*  **-  **-  ^*  V*  *•*  5&  ?fe  -J^  *JU  V*  ^  *11'  ¥*  'J'*  *'  ^*  ,*1-  ***  ^r  ■sV  ^  ^^  ^^  **-  -J-  ***  *** 
^*  T1  ^  n5  3jfc  ji*  ^i-  ^r  ^*  t-  -?  *F  *.■»  -  -  *-*■  «v  'i*  n*  ~*  *i*  ~r*  "v-  t  t*  -i*  *?•  T"  t*  -<-  n*  ^i*  o-  n*  i*  *i*  -r*  *-  -<*  t-  *c-  t*  -r-  v  i*  *.*  *i*  ^  t  *r*  'r*  -v  -»*  -i*  t- 

C  *  FORTRAN    PROGRAM    FOR    THIRD    ORDER    ADAMS    PREDICTOR-  * 

C  *  CORRECTOR    METHOD    IN    THE    NUMERICAL    SOLUTION    OF    THE  * 

C  *  ORDINARY    DIFFERENTIAL    EQUATION    IN    THE    FORM    OF  * 

C  *  DY/DX=F(X,Y)    WITH    THE     INITIAL    CONDITION    Y(X0)=Y0  * 

C  *  THE    PARAMETERS    TO    THE    PROGRAM    HAVE    THE    FOLLOWING  * 

C  *  MEANINGS...  * 

C  *         XMAX THE    TERMINAL    BOUNDARY    CONDITION  * 

C  *         EPS THE    CONVERGENCE    TEST    CONSTANT  * 

C  *         MAX THE    MAXIMUM    NUMBER    OF    ITERATIONS  * 

C  *         YP THE    PREDICTED    VALUE    OF    THE    NEXT  * 

C  *  SOLUTION    POINT  * 

C  *         YC THE    CORRECTED    VALUE    OF    THE    NEXT  * 

C  *  SOLUTION    POINT  * 

C  *         FP THE    FUNCTION    EVALUATION  * 

C  *         H__ THE    STEP    SIZE    TO    BE    USED    TO    ADVANCE  * 

C  *  THE    POINT    OF     INTEGRATION  * 

C  *         ICON THE    DESIRED    POINT    OF    PRINTING    THE  * 

C  *  COMPUTED    SOLUTION    FOP    EVERY    FIXED  * 

C  *  NUMBER    OF    INTEGRATION    POINTS  * 

C  *         YEXACT THE    TRUE    SOLUTION    VALUES  * 

C  *          ERROR m THE    ABSOLUTE    DEVIATION    OF    THE    NUMERICAL* 

C  *  SOLUTION    FROVt    THE    TRUE    SOLUTION  * 

C  *         FCT FUNCTION    SUBROUTINE    TO    EVALUATE    ODE  * 

C  *  FUNCTION  * 

C  *         EXACT FUNCTION    SUBROUTINE    TO    COMPUTE    TRUE  * 

C  *  SOLUTION    VALUES  * 

C  *         RKUTTA THE    SUBROUTINE    USED    TO    GENERATE  * 

C  *  NEEDED     STARTING    VALUES  * 

C 

C 

C  INITIALIZE    INPUT    CONSTANTS 

C 

XMAX=10.0 

MAX=2 

EPS=0. OOOOOl 

K=2 
C 

C  READ    STEP    SIZE    AND    CONTROL    PRINTING 

C 

5  READ(5,100)    H,ICON 

C 
C  SPECIFY    INITIAL    CONDITION 

X0=0.0 

Y0=1.0 

YR  =  YO 
C 

C  TEST    FOR    END    OF    INPUT    DATA 

C 

IF(ICON-O)    45,45,7 
C 

C  PRINT    INPUT    DATA    SET 

C 
7  WRITE(6,1000)    H, ICON 

WRITE(6,1001)    XO,YO 
C 

C  COMPUTE    NEXT    SOLUTION    POINT 

C 
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NC=( XMAX+H/2.J/H 

FO=FCT(XO,YO) 

X1=X0+H 

X2=X0+2.*H 

X=X0+3.*H 

CALL    RKUTTA(X0,YR,Y1,Y2,K,H) 

F1=FCT(X1,Y1) 

F2=FCT(X2,Y2) 

DO    40    N=3,NC 

M=l 

YP=Y1+H/12.*(23.*F2-16.*F1+5.*F0) 

FP=FCT(X,YP) 
11  YC  =  Y2+H/12.*(5.*FP+8.*F2-tFl) 

DELY=ABS(YC-YP) 
C 

C  TEST    FOR    CONVERGENCE 

C 

IF(DELY-EPS)    30,30,15 
C 

C      TEST  IF  MAXIMUM  NUMBER  OF  ITERATIONS  IS  SATISFIED 
C 
15    IF(M-MAX)  20,20,30 
20    FP=FCT(X,YC) 

M=M+1 

GO  TO  11 
C 

C  COMPUTE    TRUE    SOLUTION    VALUE 

C 

30         YEXACT=EXACT(X) 
C 

C  COMPUTE    ERROR    IN    CALCULATION 

C 

ERROR=YEXACT-YC 
C 

C  TEST    IF    DESIRED    POINT    OF    PRINTING    IS    REACHED 

C 

IF( (N/ICON*ICON).EQ.N)    WRIT E (6 , 1002)    X, YP , YC , YEXACT, 
1ERR0R 
C 

C  UPDATE    REQUIRED    PARAMETER    VALUES 

C 

F0=F1 

F1=F2 

F2=FP 

Y0=Y1 

Y1=Y2 

Y2=YC 

X=X+H 
40         CONTINUE 
C 

C  LOOP    BACK    TO    READ    NEXT    SET    OF    INPUT    DATA 

C 

GO    TO    5 
45         STOP 
100       F0RMAT(FL0.7,I10) 

1000  FORMAT ( ////4 IX , ' *  I NPUT    DATA    SET    US  ED* • //43X,  ■  H  =' 
1F10.7/43X, "ICON    =■ ,1 10) 

1001  FORMAT  (//29X,«**RESULTS    OF    ADAM3    PREDICTOR-CORRECTOR 
1METH0D**I//10X,»  X       ■ f5X,«  YP  •  ,5X,'  YC  ■ 

15X,'  YEXACT  «,5X,«  ERROR  « //12 X,F1 0.7 ,44X , 

1F10.7) 

1002  FORMAT ( '0,,9X,F15.7,2X,F15.7,2X,F15.7,2X,F15.7,2X, 
1F15.7) 

END 

SUBROUTINE    RKUTTA( XX , YY, Yl , Y2 , KK,HH) 
NN  =  1 
17  CK1=HH*FCT(XX,YY) 

CK2=HH*FCT(XX+HH/2.0,YY+CK 1/2.0) 
CK3=HH*FCT(XX+HH/2.0,YY+CK2/2.0) 
CK4-HH*FCT(XX+HHfYY+CK3) 
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GO  TO  {  101,102) ,NN 

101  Yl=YY+(CKl+2.0*CK2+2.0*CK3+CK4)/6.0 
YY=Y1 

77    XX=XX+HH 

YEXACT=EXACT(XX) 

ERROR=YEXACT-YY 

WRITE(6,1004)  XX, YY,YEXACT, ERROR 

NN=NN+1 

IF{NN-KK)  17,17,99 

102  Y2=YY+(CKl+2.0*CK2+2.0*CK3+CK4)/6.0 
YY  =  Y2 

GO  TO  77 
1004  FORMAT! ,0',9X,F15.7,2X,F15.7,16X,F15.7,2X,F15.7) 
99    RETURN 

END 

FUNCTION  FCT(XN,YN) 

FCT=-YN 

RETURN 

END 

FUNCTION  EXACT(XK) 

EXACT=EXP(-XK) 

RETURN 

END 
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*INPUT  DATA  SET  USED* 
H     =  0.5000000 
ICON  =  2 

**RESULTS  OF  ADAM3  PREDICTOR-CORRECTOR  METHOD** 


0.0000000 
0.5000000 
1.0000000 
2.0000000 
3.0000000 
4.0000000 
6.0000000 
7.0000000 
8.0000000 
9.0000000 
10.0000000 


YP 

0.60677  09 
0.3681710 
0.2630881 
0.0949250 
0.0333163 
0.0041118 
0.0014410 
0.0005048 
0.0001768 
0.0000619 


YC 


0.1307259 
0.0457278 
0.0160214 
0.0019661 
0.0006887 
0.0002412 
0.0000845 
0.0000296 


YEXACT 
1.0000000 
0.6065306 
0.3678794 
0.1353353 
0.0497871 
0.0183156 
0.0024788 
0.0009119 
0.0003355 
0.0001234 
0.0000454 


ERROR 

-0.0002403 
-0.0002916 
0.0046093 
0.0040593 
0.0022942 
0.0005126 
0.0002232 
0.0000942 
0.0000389 
0.0000158 


*INPUT  DATA  SET  USED* 
H     =  1.0000000 
ICON  =  1 


**RESULTS  OF  ADAM3 

PREDICTOR- 

•CORRECTOR  METHOD** 

X 

YP 

YC 

YEXACT 

ERROR 

0, 

,0000000 

1.0000000 

1. 

,0000000 

0.3750000 

0.3678794 

-0.0071206 

2. 

,0000000 

0.1406250 

0.1353353 

-0.0052897 

3. 

,0000000 

0.1888022 

0.0454788 

0.0497871 

0.0043082 

4, 

,0000000 

0.0217022 

0.0021872 

0.0183156 

0.0161284 

5. 

,0000000 

0.0785821 

-0.0024490 

0.0067379 

0.0091869 

6. 

,0000000 

-0.0525024 

-0.0057783 

0.0024788 

0.0082570 

7. 

,0000000 

-0.0479046 

0.0015025 

0.0009119 

-0.0005906 

8, 

,0000000 

-0.0577729 

-0.0018528 

0.0003355 

0.0021883 

9. 

,0000000 

0.0527026 

0.0029584 

0.0001234 

-0.0028350 

0, 

,0000000 

-0.0540226 

-0.0020290 

0.00004  54 

0.0020744 
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PROGRAM  EIGHT 

FOURTH  ORDER  ADAMS  PREDICTOR-CORRECTOR  METHOD.  ODE 
Y«=-Y  WITH  Y(0)=1.  TO  SOLVE  ANOTHER  ODE  CHANGE 
INITIAL  CONDITIONS  AND  FUNCTION  SUBROUTINES. 


C  *  FORTRAN  PROGRAM  FOR  FOURTH  ORDER  ADAMS  PREDICTOR-  * 

C  *  CORRECTOR  METHOD  IN  THE  NUMERICAL  SOLUTION  OF  THE  * 

C  *  ORDINARY  DIFFERENTIAL  EQUATION  IN  THE  FORM  OF  * 

C  *  DY/DX=F(X,Y)  WITH  THE  INITIAL  CONDITION  Y(XO)=YO  * 

C  *  THE  PARAMETERS  TO  THE  PRCGRAM  HAVE  THE  FGLLOWING  * 

C  *  MEANINGS...  * 

C  *    XMAX THE  TERMINAL  BOUNDARY  CONDITION  * 

C  *    EPS THE  CONVERGENCE  TEST  CONSTANT  * 

C  *    MAX THE  MAXIMUM  NUMBER  OF  ITERATIONS  * 

C  *    YP "THE  PREDICTED  VALUE  OF  THE  NEXT  * 

C  *  SOLUTION  POINT  * 

C  *    YC THE  CORRECTED  VALUE  OF  THE  NEXT  * 

C  *  SOLUTION  POINT  * 

C  *    FP THE  FUNCTION  EVALUATION  * 

C  *    H THE  STEP  SIZE  TO  BE  USED  TO  ADVANCE  * 

C  *  ""THE  POINT  OF  INTEGRATION  * 

C  *    ICON THE  DESIRED  POINT  OF  PRINTING  THE  * 

C  *  COMPUTED  SOLUTION  FOR  EVERY  FIXED  * 

C  *  NUMBER  OF  INTEGRATION  POINTS  * 

C  *    YEXACT THE  TRUE  SOLUTION  VALUES  * 

C  *    ERROR I_THE  ABSOLUTE  DEVIATION  OF  THE  NUMERICAL* 

C  *  SOLUTION  FROM  THE  TRUE  SOLUTION  * 

C  *    FCT  FUNCTION  SUBROUTINE  TO  EVALUATE  ODE  * 

C  *  FUNCTION  * 

C  *    EXACT FUNCTION  SUBROUTINE  TO  COMPUTE  TRUE  * 

C  *  SOLUTION  VALUES  * 

C  *    RKUTTA THE  SUBROUTINE  USED  TO  GENERATE  * 

C  *^  ^^    DEEDED  STARTING ^VALUES^  * 

C 

c 

C      INITIALIZE  INPUT  CONSTANTS 
C 

XMAX=10.0 

MAX  =  2 

EPS=0. OO0001 

K=3 
C 

C      READ  STEP  SIZE  AND  CONTROL  PRINTING 
C 

5     READ(5,100)  H,ICON 
C 

C  SPECIFY    INITIAL    CONDITION 

C 

X0=0.0 

Y0=1.0 

YR=Y0 
C 

C      TEST  FOR  END  OF  INPUT  DATA 
C 

IF(ICON-O)  45,45,7 

C      PRINT  INPUT  DATA  SET 
C 
7     WRITE(6,1000)  H,ICON 
WRITE(6, 1001)  X0,Y0 
C 
C      COMPUTE  NEXT  SOLUTION  POINT 
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NC=(XMAX+H/2.)/H 

F0=FCT(XO,Y0) 

X1=X0+H 

X2=X0+2.*H 

X3=X0+3.*H 

X=X0+4.*H 

CALL    RKUTTA(X0,YR,Y1,Y2,Y3,K,H) 

Fl=FCT(XltYl) 

F2=FCT(X2,Y2) 

F3=FCT(X3,Y3) 

DO    40    N=4,NC 

M=l 

YP=Y3+H/24.*(5  5.*F3-59.*F2+37.*F1-9.*F0) 

FP=FCT(X,YP) 
11         YC=Y3+H/24.*(9.*FP+19.*F3-5.*F2+F1) 

DELY=ABS(YC-YP) 
C 

C      TEST  FOR  CONVERGENCE 
C 

IF(DELY-EPS)  30,30,15 
C 

C  TEST     IF    MAXIMUM    NUMBER    OF    ITERATIONS    IS    SATISFIED 

C 
15  IF(M-MAX)     20,20,30 

20  FP=FCT(X,YC) 

M=M+1 

GO    TO    11 
C 

C  COMPUTE    TRUE    SOLUTION    VALUE 

C 

30         YEXACT=EXACT(X) 
C 

C  COMPUTE    ERROR    IN    CALCULATION 

C 

ERROR=YEXACT-YC 
C 

C  TEST    IF    DESIRED    POINT    OF    PRINTING    IS    REACHED 

C 

IF( (N/ICON*ICON).EQ.N)     WRITE (6 , 1 002 )    X , YP, YC, YEXACT , 
1ERR0R 
C 

C  UPDATE    REQUIRED    PARAMETER    VALUES 

C 

Y0=Y1 

Y1=Y2 

Y2=Y3 

Y3=YC 

F0  =  F1 

F1  =  F2 

F2=F3 

F3=FP 

40         CONTINUE 
C 

C  LOOP    BACK    TO    READ    NEXT    SET    OF    INPUT    DATA 

C 

GO    TO    5 

45         STOP 

100       FORMAT(F10.7,I10) 

1000  F0RMATI////41X, «*INPUT  DATA  SET  USED*' //43X , 'H     =  • 
1F10.7/43X," ICON  =•  ,110) 

1001  FORMAT (//29X,»**RESULTS  OF  ADAM4  PREDICTOR-CORRECTOR 
1METHCD**,//10X, •    X   ',5X,»    YP    t,5X11  YC    • 

15X,'    YEXACT    '^X,*    ERROR    •  //  12X,  F10  .7  ,44X  , 
1F10.7) 

1002  FORMAT  CO'  ,9X,  F15.7  ,  2X  ,  F  15  .7  ,  2X,  F15.  7,  2X,  F  15.  7,  2X, 
1F15.7) 

END 
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SUBROUTINE  RKUTTA( XX , YY , Yl , Y2 , Y3 ,KK, HH) 
NN=1 
17    CK1=HH*FCT(XX,YY) 

CK2=Hh*FCT(XX+HH/2.0,YY+CK 1/2.0) 
CK3=HH*FCT(XX+HH/2.0,YY+CK 2/2.0) 
CK4=HH*FCT(XX+HH,YY+CK3) 
GO  TO  ( 101,102,103) ,NN 

101  Yl=YY+(CKl+2.0*CK2+2.0*CK3+CK4)/6.0 
YY=Y1 

77    XX=XX+HH 

YEXACT=EXACT(XX) 

ERROR=YEXACT-YY 

WRITE(6,1004)  XX, YY,YEXACT, ERROR 

NN=NN+1 

IF(NN-KK)  17,17,99 

102  Y2=YY+(CKl+2.0*CK2+2.0*CK3+CK4)/6.0 
YY=Y2 

GO  TO  77 

103  Y3=YY+(CKl+2.0*CK2+2.0*CK3+CK4)/6. 0 
YY=Y3 

GO  TO  77 
1004  FORMAT (  ,0',9X,F15.7,2X,F15.7,16X,F15.7,2X,F15.7) 
99    RETURN 

END 

FUNCTION  FCT(XN,YN) 

FCT=-YN 

RETURN 

END 

FUNCTION  EXACT(XK) 
EXACT=EXP(-XK) 

RETURN 
END 
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*INPUT  DATA  SET  USED* 
H     =  0.5000000 
ICON  =  2 

**RESULTS  OF  ADAM4  PREDICTOR-CORRECTOR  METHOD** 


0.0000000 
0.5000000 
1.0000000 
1.5000000 
2.0000000 
3.0000000 
4.0000000 
5.0000000 
6.0000000 
7.0000000 
8.0000000 
9.0000000 
10.0000000 


YP 

0.6067709 
0.3681710 
0.2233955 
0.1397457 
0.0512678 
0.0187947 
0.0068879 
0.0025232 
0.0009245 
0.0003387 
0.0001241 
0.0000455 


YC 


0.1352788 
0.0495746 
0.0181669 
0.0066569 
0.0024393 
0.0008938 
0.0003275 
0.0001200 
0.0000440 


YEXACT 
1.0000000 
0.6065306 
0.3678794 
0.2231301 
0.1353353 
0.0497871 
0.0183156 
0.0067379 
0.0024788 
0.0009119 
0.0003355 
0.0001234 
0.0000454 


ERROR 

-0.0002403 
-0.0002916 
-0.0002654 
0.0000565 
0.0002124 
0.0001488 
0.0000810 
0.0000394 
0.0000180 
0.0000079 
0.0000034 
0.0000014 


*INPUT  DATA  SET  USED* 
H     =  1.0000000 
ICON  =  1 

**RESULTS  OF  ADAM4  PREDICTOR-CORRECTOR  METHOD** 


0.0000000 
1.0000000 
2.0000000 
3.0000000 
4.0000000 
5.0000000 
6.0000000 
7.0000000 
8.0000000 
9.0000000 
10.0000000 


YP 

0.3750000 
C. 1406250 
0.0527344 
0.0744628 
0.0091045 
0.0319238 

-0.0306429 
0.0368799 

-0.0442679 
0.0550162 


YC 


0.0149522 
-0.0007950 
-0.00Q4662 
-0.0027268 

0.0015700 
-0.0027738 

0.0028112 


YEXACT 
1.0000000 
0.3678794 
0.1353353 
0.0497871 
0.0183156 
0.0067379 
0.0024788 
0.0009119 
0.0003355 
0.0001234 
0.0000454 


ERROR 

-0.0071206 

-0.0052897 

-0.0029473 

0.0033634 

0.0075330 

0.0029450 

0.0036387 

-0.0012345 

0.0028972 

-0.0027658 
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